use std::sync::Arc;
use antecedent_core::PriorAssumption;
use crate::error::ProbError;
use crate::prior::{GaussianCoefficientPrior, PriorSet, PriorSpec};
const COMPOSE_VAR_FLOOR: f64 = 1e-12;
#[must_use]
pub fn kish_ess(weights: &[f64]) -> f64 {
let sum: f64 = weights.iter().sum();
let sum_sq: f64 = weights.iter().map(|w| w * w).sum();
if sum_sq > 0.0 { (sum * sum) / sum_sq } else { 0.0 }
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct ExternalPriorWeight {
pub alpha: f64,
pub mixture_weight: Option<f64>,
}
impl ExternalPriorWeight {
pub fn new(alpha: f64, mixture_weight: Option<f64>) -> Result<Self, ProbError> {
let w = Self { alpha, mixture_weight };
w.validate()?;
Ok(w)
}
pub fn power(alpha: f64) -> Result<Self, ProbError> {
Self::new(alpha, None)
}
pub fn power_mixture(alpha: f64, mixture_weight: f64) -> Result<Self, ProbError> {
Self::new(alpha, Some(mixture_weight))
}
pub fn validate(self) -> Result<(), ProbError> {
if !self.alpha.is_finite() || !(0.0..=1.0).contains(&self.alpha) {
return Err(ProbError::InvalidPrior {
message: "external prior alpha must be finite and in [0, 1]",
});
}
if let Some(w) = self.mixture_weight {
if !w.is_finite() || w < 0.0 {
return Err(ProbError::InvalidPrior {
message: "mixture weight must be finite and >= 0",
});
}
}
Ok(())
}
}
#[derive(Clone, Debug, PartialEq)]
pub struct ExternalPriorSource {
pub id: Arc<str>,
pub prior: PriorSet,
pub weight: ExternalPriorWeight,
pub ess: Option<f64>,
}
impl ExternalPriorSource {
pub fn validate(&self) -> Result<(), ProbError> {
self.weight.validate()?;
if let Some(ess) = self.ess {
if !ess.is_finite() || ess < 0.0 {
return Err(ProbError::InvalidPrior {
message: "external prior source ess must be finite and >= 0",
});
}
}
Ok(())
}
}
#[derive(Clone, Debug, PartialEq)]
pub struct ComposedPrior {
pub prior: PriorSet,
pub source_ids: Arc<[Arc<str>]>,
pub alphas_requested: Arc<[f64]>,
pub alphas_applied: Arc<[f64]>,
pub mixture_weights: Arc<[Option<f64>]>,
pub effective_ess: Arc<[Option<f64>]>,
pub composed_ess: Option<f64>,
pub kish_ess: Option<f64>,
}
impl ComposedPrior {
#[must_use]
pub fn as_prior_set(&self) -> &PriorSet {
&self.prior
}
#[must_use]
pub fn into_prior_set(self) -> PriorSet {
self.prior
}
}
pub fn compose_external_priors(
sources: &[ExternalPriorSource],
baseline: &PriorSet,
) -> Result<ComposedPrior, ProbError> {
let alphas: Vec<f64> = sources.iter().map(|s| s.weight.alpha).collect();
compose_external_priors_with_alphas(sources, &alphas, &alphas, baseline)
}
pub fn compose_external_priors_with_alphas(
sources: &[ExternalPriorSource],
alphas_requested: &[f64],
alphas_applied: &[f64],
baseline: &PriorSet,
) -> Result<ComposedPrior, ProbError> {
if sources.len() != alphas_requested.len() || sources.len() != alphas_applied.len() {
return Err(ProbError::Shape {
message: "compose_external_priors: alpha vector length mismatch",
});
}
for &a in alphas_requested.iter().chain(alphas_applied.iter()) {
if !a.is_finite() || !(0.0..=1.0).contains(&a) {
return Err(ProbError::InvalidPrior {
message: "external prior alpha must be finite and in [0, 1]",
});
}
}
for src in sources {
src.validate()?;
}
validate_mixture_weights(sources)?;
let base_coef = baseline.gaussian_coefficients().ok_or(ProbError::InvalidPrior {
message: "baseline prior missing GaussianCoefficients",
})?;
base_coef.validate()?;
let n = base_coef.len();
for src in sources {
let coef = src.prior.gaussian_coefficients().ok_or(ProbError::InvalidPrior {
message: "external source prior missing GaussianCoefficients",
})?;
coef.validate()?;
if coef.len() != n {
return Err(ProbError::Shape {
message: "compose_external_priors: coefficient dimension mismatch",
});
}
}
let use_mixture = sources.iter().any(|s| s.weight.mixture_weight.is_some());
let composed_coef = if use_mixture {
compose_mixture(base_coef, sources, alphas_applied)?
} else {
compose_power_add(base_coef, sources, alphas_applied)?
};
let mut prior = PriorSet {
specs: Vec::new(),
contrast: baseline.contrast,
categorical: baseline.categorical.clone(),
restrictions: Vec::new(),
};
prior.push(PriorSpec::GaussianCoefficients(composed_coef));
if let Some(ig) = baseline.residual_inv_gamma() {
prior.push(PriorSpec::ResidualInvGamma(ig));
} else if let Some(v) = baseline.known_residual_variance() {
prior.push(PriorSpec::KnownResidualVariance(v));
}
for r in &baseline.restrictions {
prior.restrictions.push(r.clone());
}
for src in sources {
for r in &src.prior.restrictions {
prior.restrictions.push(r.clone());
}
}
prior.restrictions.push(composition_assumption(sources, alphas_requested, alphas_applied));
prior.validate()?;
let source_ids: Vec<Arc<str>> = sources.iter().map(|s| Arc::clone(&s.id)).collect();
let mixture_weights: Vec<Option<f64>> =
sources.iter().map(|s| s.weight.mixture_weight).collect();
let effective_ess = effective_ess_per_source(sources, alphas_applied, use_mixture);
let composed_ess = if use_mixture { None } else { power_composed_ess(sources, alphas_applied) };
let kish_ess_diag = if sources.is_empty() {
None
} else {
Some(kish_ess(&kish_weights_for_composition(sources, alphas_applied, use_mixture)))
};
Ok(ComposedPrior {
prior,
source_ids: Arc::from(source_ids),
alphas_requested: Arc::from(alphas_requested.to_vec()),
alphas_applied: Arc::from(alphas_applied.to_vec()),
mixture_weights: Arc::from(mixture_weights),
effective_ess: Arc::from(effective_ess),
composed_ess,
kish_ess: kish_ess_diag,
})
}
fn effective_ess_per_source(
sources: &[ExternalPriorSource],
alphas_applied: &[f64],
use_mixture: bool,
) -> Vec<Option<f64>> {
sources
.iter()
.zip(alphas_applied.iter())
.map(|(src, &alpha)| {
let dropped = if use_mixture {
alpha <= 0.0 || src.weight.mixture_weight.unwrap_or(0.0) <= 0.0
} else {
alpha == 0.0
};
src.ess.map(|ess| if dropped { 0.0 } else { alpha * ess })
})
.collect()
}
fn power_composed_ess(sources: &[ExternalPriorSource], alphas_applied: &[f64]) -> Option<f64> {
let mut total = 0.0;
for (src, &alpha) in sources.iter().zip(alphas_applied.iter()) {
if alpha == 0.0 {
continue;
}
total += alpha * src.ess?;
}
Some(total)
}
fn kish_weights_for_composition(
sources: &[ExternalPriorSource],
alphas_applied: &[f64],
use_mixture: bool,
) -> Vec<f64> {
if use_mixture {
sources
.iter()
.zip(alphas_applied.iter())
.map(
|(src, &alpha)| {
if alpha <= 0.0 { 0.0 } else { src.weight.mixture_weight.unwrap_or(0.0) }
},
)
.collect()
} else {
alphas_applied.to_vec()
}
}
fn validate_mixture_weights(sources: &[ExternalPriorSource]) -> Result<(), ProbError> {
if sources.is_empty() {
return Ok(());
}
let any = sources.iter().any(|s| s.weight.mixture_weight.is_some());
let all = sources.iter().all(|s| s.weight.mixture_weight.is_some());
if any && !all {
return Err(ProbError::InvalidPrior {
message: "mixture weights must be set on all sources or none",
});
}
if !any {
return Ok(());
}
let sum: f64 = sources.iter().map(|s| s.weight.mixture_weight.unwrap_or(0.0)).sum();
if !sum.is_finite() || sum > 1.0 + 1e-12 {
return Err(ProbError::InvalidPrior { message: "sum of mixture weights must be <= 1" });
}
Ok(())
}
fn compose_power_add(
baseline: &GaussianCoefficientPrior,
sources: &[ExternalPriorSource],
alphas: &[f64],
) -> Result<GaussianCoefficientPrior, ProbError> {
let n = baseline.len();
let mut lam = baseline.precision();
let mut num = vec![0.0; n];
for i in 0..n {
num[i] = lam[i] * baseline.mean[i];
}
for (src, &alpha) in sources.iter().zip(alphas.iter()) {
if alpha == 0.0 {
continue;
}
let coef = src.prior.gaussian_coefficients().expect("validated");
let prec = coef.precision();
for i in 0..n {
let a_lam = alpha * prec[i];
lam[i] += a_lam;
num[i] += a_lam * coef.mean[i];
}
}
let mut mean = vec![0.0; n];
let mut variance = vec![0.0; n];
for i in 0..n {
if !(lam[i] > 0.0) || !lam[i].is_finite() {
return Err(ProbError::Numerical {
message: "compose_external_priors: non-positive composed precision".into(),
});
}
mean[i] = num[i] / lam[i];
variance[i] = (1.0 / lam[i]).max(COMPOSE_VAR_FLOOR);
}
let out = GaussianCoefficientPrior { mean: Arc::from(mean), variance: Arc::from(variance) };
out.validate()?;
Ok(out)
}
fn compose_mixture(
baseline: &GaussianCoefficientPrior,
sources: &[ExternalPriorSource],
alphas: &[f64],
) -> Result<GaussianCoefficientPrior, ProbError> {
let n = baseline.len();
let mut active_w = 0.0;
let mut comps: Vec<(f64, &GaussianCoefficientPrior, f64)> = Vec::new();
for (src, &alpha) in sources.iter().zip(alphas.iter()) {
let w = src.weight.mixture_weight.unwrap_or(0.0);
if alpha <= 0.0 || w <= 0.0 {
continue;
}
let coef = src.prior.gaussian_coefficients().expect("validated");
comps.push((w, coef, alpha));
active_w += w;
}
let leftover = (1.0 - active_w).max(0.0);
if leftover > 0.0 {
comps.push((leftover, baseline, 1.0));
}
if comps.is_empty() {
return Err(ProbError::InvalidPrior {
message: "compose_external_priors: mixture has no positive-mass components",
});
}
let mut mean = vec![0.0; n];
let mut variance = vec![0.0; n];
for i in 0..n {
let mut mu = 0.0;
let mut second = 0.0;
for &(w, coef, alpha) in &comps {
let m = coef.mean[i];
let v = (coef.variance[i] / alpha).max(COMPOSE_VAR_FLOOR);
mu += w * m;
second += w * (v + m * m);
}
mean[i] = mu;
variance[i] = (second - mu * mu).max(COMPOSE_VAR_FLOOR);
}
let out = GaussianCoefficientPrior { mean: Arc::from(mean), variance: Arc::from(variance) };
out.validate()?;
Ok(out)
}
fn composition_assumption(
sources: &[ExternalPriorSource],
alphas_requested: &[f64],
alphas_applied: &[f64],
) -> PriorAssumption {
let mut parts = Vec::with_capacity(sources.len());
for (i, src) in sources.iter().enumerate() {
let w = src.weight.mixture_weight.map_or_else(|| "none".to_string(), |x| format!("{x}"));
parts.push(format!(
"{}:alpha_req={},alpha_app={},w={}",
src.id, alphas_requested[i], alphas_applied[i], w
));
}
PriorAssumption {
id: Arc::from("external_composed_prior"),
description: Arc::from(format!(
"External power-prior / mixture composition [{}]",
parts.join("; ")
)),
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::prior::GaussianCoefficientPrior;
fn gauss(mean: f64, var: f64) -> PriorSet {
let mut p = PriorSet::new();
p.push(PriorSpec::GaussianCoefficients(
GaussianCoefficientPrior::shared(1, mean, var).unwrap(),
));
p
}
#[test]
fn rejects_alpha_out_of_range() {
assert!(ExternalPriorWeight::power(-0.1).is_err());
assert!(ExternalPriorWeight::power(1.1).is_err());
assert!(ExternalPriorWeight::power(f64::NAN).is_err());
}
#[test]
fn rejects_mixture_weight_sum_gt_one() {
let baseline = PriorSet::weakly_informative(1);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(1.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.7).unwrap(),
ess: None,
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.5).unwrap(),
ess: None,
},
];
let err = compose_external_priors(&sources, &baseline).unwrap_err();
assert!(matches!(err, ProbError::InvalidPrior { .. }));
}
#[test]
fn power_prior_precision_add_analytic() {
let baseline = gauss(0.0, 4.0);
let sources = [ExternalPriorSource {
id: Arc::from("old"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(0.5).unwrap(),
ess: None,
}];
let composed = compose_external_priors(&sources, &baseline).unwrap();
let coef = composed.prior.gaussian_coefficients().unwrap();
let lam = 1.0 / coef.variance[0];
assert!((lam - 0.75).abs() < 1e-12);
assert!((coef.mean[0] - (4.0 / 3.0)).abs() < 1e-12);
assert!(composed.prior.restrictions.iter().any(|r| &*r.id == "external_composed_prior"));
}
#[test]
fn mixture_preserves_leftover_baseline_mass() {
let baseline = gauss(0.0, 100.0);
let sources = [ExternalPriorSource {
id: Arc::from("s"),
prior: gauss(10.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.4).unwrap(),
ess: None,
}];
let composed = compose_external_priors(&sources, &baseline).unwrap();
let coef = composed.prior.gaussian_coefficients().unwrap();
assert!((coef.mean[0] - 4.0).abs() < 1e-10);
assert!((coef.variance[0] - 84.4).abs() < 1e-10);
}
#[test]
fn applied_alpha_override() {
let baseline = gauss(0.0, 4.0);
let sources = [ExternalPriorSource {
id: Arc::from("old"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(1.0).unwrap(),
ess: None,
}];
let composed =
compose_external_priors_with_alphas(&sources, &[1.0], &[0.0], &baseline).unwrap();
let coef = composed.prior.gaussian_coefficients().unwrap();
assert!((coef.mean[0] - 0.0).abs() < 1e-12);
assert!((coef.variance[0] - 4.0).abs() < 1e-12);
assert_eq!(&*composed.alphas_requested, &[1.0]);
assert_eq!(&*composed.alphas_applied, &[0.0]);
}
#[test]
fn rejects_mixed_mixture_mode() {
let baseline = PriorSet::weakly_informative(1);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(1.0, 1.0),
weight: ExternalPriorWeight::power(1.0).unwrap(),
ess: None,
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.5).unwrap(),
ess: None,
},
];
assert!(compose_external_priors(&sources, &baseline).is_err());
}
#[test]
fn power_prior_ess_sums_over_contributing_sources() {
let baseline = gauss(0.0, 4.0);
let sources = [ExternalPriorSource {
id: Arc::from("old"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(0.5).unwrap(),
ess: Some(40.0),
}];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert_eq!(composed.effective_ess.len(), 1);
assert!((composed.effective_ess[0].unwrap() - 20.0).abs() < 1e-12);
assert!((composed.composed_ess.unwrap() - 20.0).abs() < 1e-12);
assert!((composed.kish_ess.unwrap() - 1.0).abs() < 1e-12);
}
#[test]
fn power_prior_composed_ess_none_without_full_ess_coverage() {
let baseline = gauss(0.0, 4.0);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(0.5).unwrap(),
ess: Some(40.0),
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(3.0, 1.0),
weight: ExternalPriorWeight::power(0.25).unwrap(),
ess: None,
},
];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert!((composed.effective_ess[0].unwrap() - 20.0).abs() < 1e-12);
assert!(composed.effective_ess[1].is_none());
assert!(composed.composed_ess.is_none());
assert!((composed.kish_ess.unwrap() - 1.8).abs() < 1e-12);
}
#[test]
fn power_prior_dropped_source_contributes_no_ess() {
let baseline = gauss(0.0, 4.0);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(0.5).unwrap(),
ess: Some(40.0),
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(5.0, 1.0),
weight: ExternalPriorWeight::power(0.0).unwrap(),
ess: Some(999.0),
},
];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert!((composed.effective_ess[0].unwrap() - 20.0).abs() < 1e-12);
assert!((composed.effective_ess[1].unwrap() - 0.0).abs() < 1e-12);
assert!((composed.composed_ess.unwrap() - 20.0).abs() < 1e-12);
assert!((composed.kish_ess.unwrap() - 1.0).abs() < 1e-12);
}
#[test]
fn mixture_composed_ess_is_none_but_effective_ess_reported_per_source() {
let baseline = gauss(0.0, 100.0);
let sources = [ExternalPriorSource {
id: Arc::from("s"),
prior: gauss(10.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.4).unwrap(),
ess: Some(50.0),
}];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert!(composed.composed_ess.is_none());
assert!((composed.effective_ess[0].unwrap() - 50.0).abs() < 1e-12);
assert!((composed.kish_ess.unwrap() - 1.0).abs() < 1e-12);
}
#[test]
fn mixture_dropped_source_contributes_no_effective_ess() {
let baseline = gauss(0.0, 100.0);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(10.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.4).unwrap(),
ess: Some(50.0),
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(20.0, 1.0),
weight: ExternalPriorWeight::power_mixture(1.0, 0.0).unwrap(),
ess: Some(999.0),
},
];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert!((composed.effective_ess[0].unwrap() - 50.0).abs() < 1e-12);
assert!((composed.effective_ess[1].unwrap() - 0.0).abs() < 1e-12);
assert!(composed.composed_ess.is_none());
assert!((composed.kish_ess.unwrap() - 1.0).abs() < 1e-12);
}
#[test]
fn sources_without_ess_report_none_effective_and_composed() {
let baseline = gauss(0.0, 4.0);
let sources = [
ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(2.0, 1.0),
weight: ExternalPriorWeight::power(0.5).unwrap(),
ess: None,
},
ExternalPriorSource {
id: Arc::from("b"),
prior: gauss(3.0, 1.0),
weight: ExternalPriorWeight::power(0.3).unwrap(),
ess: None,
},
];
let composed = compose_external_priors(&sources, &baseline).unwrap();
assert!(composed.effective_ess.iter().all(Option::is_none));
assert!(composed.composed_ess.is_none());
let kish = composed.kish_ess.unwrap();
assert!(kish.is_finite() && kish > 0.0);
}
#[test]
fn rejects_negative_ess() {
let sources = [ExternalPriorSource {
id: Arc::from("a"),
prior: gauss(1.0, 1.0),
weight: ExternalPriorWeight::power(1.0).unwrap(),
ess: Some(-1.0),
}];
let err = sources[0].validate().unwrap_err();
assert!(matches!(err, ProbError::InvalidPrior { .. }));
let baseline = PriorSet::weakly_informative(1);
assert!(compose_external_priors(&sources, &baseline).is_err());
}
#[test]
fn kish_ess_matches_transport_adjustment_formula() {
use crate::transport::TransportAdjustment;
let adj = TransportAdjustment::new([1.0, 2.0, 3.0], [0.5, 0.25, 0.25]).unwrap();
assert!((kish_ess(&[0.5, 0.25, 0.25]) - adj.kish_ess()).abs() < 1e-12);
}
}