1use std::sync::Arc;
10
11use thiserror::Error;
12
13use antecedent_core::PriorAssumption;
14
15use crate::error::ProbError;
16use crate::external_prior::{
17 ComposedPrior, ExternalPriorSource, compose_external_priors_with_alphas,
18};
19use crate::prior::{GaussianCoefficientPrior, PriorSet, PriorSpec};
20
21pub const TRANSPORT_ASSUMPTION_ID: &str = "external_transport_prior";
23
24pub const POPULATION_TAG_KEY: &str = "population";
26
27const REWEIGHT_VAR_FLOOR: f64 = 1e-12;
29
30#[derive(Clone, Copy, Debug, Eq, PartialEq, Hash)]
35pub enum TransportPolicy {
36 InvariantConditionalOutcome,
38 InvariantEffectModifiers,
40 InvariantPropensity,
43}
44
45impl TransportPolicy {
46 #[must_use]
48 pub const fn id(self) -> &'static str {
49 match self {
50 Self::InvariantConditionalOutcome => "invariant_conditional_outcome",
51 Self::InvariantEffectModifiers => "invariant_effect_modifiers",
52 Self::InvariantPropensity => "invariant_propensity",
53 }
54 }
55
56 pub fn parse(s: &str) -> Result<Self, TransportError> {
62 match s {
63 "invariant_conditional_outcome" | "InvariantConditionalOutcome" => {
64 Ok(Self::InvariantConditionalOutcome)
65 }
66 "invariant_effect_modifiers" | "InvariantEffectModifiers" => {
67 Ok(Self::InvariantEffectModifiers)
68 }
69 "invariant_propensity" | "InvariantPropensity" => Ok(Self::InvariantPropensity),
70 other => Err(TransportError::UnknownPolicy { name: Arc::from(other) }),
71 }
72 }
73}
74
75#[derive(Clone, Debug, PartialEq)]
81pub struct TransportAdjustment {
82 pub unit_effects: Arc<[f64]>,
84 pub target_weights: Arc<[f64]>,
86}
87
88impl TransportAdjustment {
89 pub fn new(
95 unit_effects: impl Into<Arc<[f64]>>,
96 target_weights: impl Into<Arc<[f64]>>,
97 ) -> Result<Self, TransportError> {
98 let unit_effects = unit_effects.into();
99 let target_weights = target_weights.into();
100 if unit_effects.is_empty() || target_weights.is_empty() {
101 return Err(TransportError::InvalidWeights {
102 message: "transport adjustment requires non-empty effects and weights",
103 });
104 }
105 if unit_effects.len() != target_weights.len() {
106 return Err(TransportError::InvalidWeights {
107 message: "unit_effects and target_weights length mismatch",
108 });
109 }
110 let mut mass = 0.0;
111 for (&e, &w) in unit_effects.iter().zip(target_weights.iter()) {
112 if !e.is_finite() {
113 return Err(TransportError::InvalidWeights {
114 message: "unit_effects must be finite",
115 });
116 }
117 if !w.is_finite() || w < 0.0 {
118 return Err(TransportError::InvalidWeights {
119 message: "target_weights must be finite and >= 0",
120 });
121 }
122 mass += w;
123 }
124 if !(mass > 0.0) {
125 return Err(TransportError::InvalidWeights {
126 message: "target_weights must have positive total mass",
127 });
128 }
129 Ok(Self { unit_effects, target_weights })
130 }
131
132 #[must_use]
134 pub fn weighted_moments(&self) -> (f64, f64) {
135 let mass: f64 = self.target_weights.iter().sum();
136 let mean = self
137 .unit_effects
138 .iter()
139 .zip(self.target_weights.iter())
140 .map(|(&e, &w)| w * e)
141 .sum::<f64>()
142 / mass;
143 let var = self
144 .unit_effects
145 .iter()
146 .zip(self.target_weights.iter())
147 .map(|(&e, &w)| {
148 let d = e - mean;
149 w * d * d
150 })
151 .sum::<f64>()
152 / mass;
153 (mean, var.max(REWEIGHT_VAR_FLOOR))
154 }
155
156 #[must_use]
161 pub fn kish_ess(&self) -> f64 {
162 crate::external_prior::kish_ess(&self.target_weights)
163 }
164}
165
166#[derive(Clone, Debug)]
168pub struct TransportContext<'a> {
169 pub source_populations: &'a [Option<&'a str>],
171 pub target_population: Option<&'a str>,
173 pub policy: Option<TransportPolicy>,
175 pub adjustment: Option<&'a TransportAdjustment>,
177 pub coef_index: Option<usize>,
179}
180
181#[derive(Clone, Debug, PartialEq)]
183pub struct TransportOutcome {
184 pub source_id: Arc<str>,
186 pub required: bool,
188 pub alpha_override: Option<f64>,
190 pub zero_reason: Option<Arc<str>>,
192}
193
194#[derive(Clone, Debug, Eq, PartialEq, Error)]
196#[non_exhaustive]
197pub enum TransportError {
198 #[error(
200 "{code}: population mismatch source={source_population:?} \
201 target={target_population:?} requires TransportPolicy",
202 code = self.code()
203 )]
204 PolicyRequired {
205 source_population: Arc<str>,
207 target_population: Arc<str>,
209 },
210 #[error(
212 "{code}: source_populations len {n_populations} != sources len {n_sources}",
213 code = self.code()
214 )]
215 SourceCountMismatch {
216 n_sources: usize,
218 n_populations: usize,
220 },
221 #[error("{code}: {message}", code = self.code())]
223 InvalidWeights {
224 message: &'static str,
226 },
227 #[error("{code}: unknown TransportPolicy `{name}`", code = self.code())]
229 UnknownPolicy {
230 name: Arc<str>,
232 },
233 #[error(
235 "{code}: coef_index {index} out of range for n_coef={n_coef}",
236 code = self.code()
237 )]
238 CoefIndexOutOfRange {
239 index: usize,
241 n_coef: usize,
243 },
244}
245
246impl TransportError {
247 #[must_use]
249 pub const fn code(&self) -> &'static str {
250 match self {
251 Self::PolicyRequired { .. } => "transport_policy_required",
252 Self::SourceCountMismatch { .. } => "transport_source_count_mismatch",
253 Self::InvalidWeights { .. } => "transport_invalid_weights",
254 Self::UnknownPolicy { .. } => "transport_unknown_policy",
255 Self::CoefIndexOutOfRange { .. } => "transport_coef_index_out_of_range",
256 }
257 }
258}
259
260impl From<TransportError> for ProbError {
261 fn from(e: TransportError) -> Self {
262 ProbError::Numerical { message: e.to_string() }
263 }
264}
265
266#[must_use]
272pub fn populations_require_transport(
273 source_population: Option<&str>,
274 target_population: Option<&str>,
275) -> bool {
276 match (source_population, target_population) {
277 (None, None) => false,
278 (Some(a), Some(b)) => a != b,
279 (Some(_), None) | (None, Some(_)) => true,
281 }
282}
283
284fn pop_label(p: Option<&str>) -> Arc<str> {
285 Arc::from(p.unwrap_or(""))
286}
287
288fn transport_assumption(
289 policy: TransportPolicy,
290 source_id: &str,
291 source_pop: Option<&str>,
292 target_pop: Option<&str>,
293 extra: &str,
294) -> PriorAssumption {
295 PriorAssumption {
296 id: Arc::from(TRANSPORT_ASSUMPTION_ID),
297 description: Arc::from(format!(
298 "TransportPolicy {} for source={source_id} source_pop={} target_pop={}{extra}",
299 policy.id(),
300 source_pop.unwrap_or(""),
301 target_pop.unwrap_or(""),
302 )),
303 }
304}
305
306fn replace_coef_moments(
307 prior: &mut PriorSet,
308 coef_index: usize,
309 mean: f64,
310 variance: f64,
311) -> Result<(), TransportError> {
312 let Some(coef) = prior.gaussian_coefficients() else {
313 return Err(TransportError::InvalidWeights {
314 message: "source prior missing GaussianCoefficients for transport reweight",
315 });
316 };
317 if coef_index >= coef.len() {
318 return Err(TransportError::CoefIndexOutOfRange { index: coef_index, n_coef: coef.len() });
319 }
320 let mut mean_v = coef.mean.to_vec();
321 let mut var_v = coef.variance.to_vec();
322 mean_v[coef_index] = mean;
323 var_v[coef_index] = variance;
324 let new_coef = GaussianCoefficientPrior { mean: Arc::from(mean_v), variance: Arc::from(var_v) };
325 new_coef.validate().map_err(|_| TransportError::InvalidWeights {
326 message: "reweighted coefficient moments invalid",
327 })?;
328 let mut specs = Vec::with_capacity(prior.specs.len());
330 let mut replaced = false;
331 for s in &prior.specs {
332 match s {
333 PriorSpec::GaussianCoefficients(_) if !replaced => {
334 specs.push(PriorSpec::GaussianCoefficients(new_coef.clone()));
335 replaced = true;
336 }
337 other => specs.push(other.clone()),
338 }
339 }
340 if !replaced {
341 specs.push(PriorSpec::GaussianCoefficients(new_coef));
342 }
343 prior.specs = specs;
344 Ok(())
345}
346
347pub fn apply_transport(
355 sources: &[ExternalPriorSource],
356 ctx: &TransportContext<'_>,
357) -> Result<(Vec<ExternalPriorSource>, Vec<TransportOutcome>), TransportError> {
358 if ctx.source_populations.len() != sources.len() {
359 return Err(TransportError::SourceCountMismatch {
360 n_sources: sources.len(),
361 n_populations: ctx.source_populations.len(),
362 });
363 }
364
365 let any_required = sources
366 .iter()
367 .zip(ctx.source_populations.iter())
368 .any(|(_, &sp)| populations_require_transport(sp, ctx.target_population));
369 if any_required && ctx.policy.is_none() {
370 for &sp in ctx.source_populations {
372 if populations_require_transport(sp, ctx.target_population) {
373 return Err(TransportError::PolicyRequired {
374 source_population: pop_label(sp),
375 target_population: pop_label(ctx.target_population),
376 });
377 }
378 }
379 }
380
381 let policy = ctx.policy;
382 let mut out_sources = Vec::with_capacity(sources.len());
383 let mut outcomes = Vec::with_capacity(sources.len());
384
385 for (src, &sp) in sources.iter().zip(ctx.source_populations.iter()) {
386 let required = populations_require_transport(sp, ctx.target_population);
387 let mut prepared = src.clone();
388 let mut alpha_override = None;
389 let mut zero_reason = None;
390
391 if required {
392 let policy = policy.expect("gated above");
393 let mut extra = String::new();
394
395 match (policy, ctx.adjustment) {
396 (TransportPolicy::InvariantPropensity, None) => {
397 alpha_override = Some(0.0);
398 zero_reason = Some(Arc::from(
399 "invariant_propensity requires target_weights; alpha forced to 0",
400 ));
401 extra.push_str("; alpha_forced=0 reason=missing_propensity_weights");
402 }
403 (_, Some(adj)) => {
404 let (mean, var) = adj.weighted_moments();
405 let n_coef = prepared
406 .prior
407 .gaussian_coefficients()
408 .ok_or(TransportError::InvalidWeights {
409 message: "source prior missing GaussianCoefficients for transport reweight",
410 })?
411 .len();
412 let idx = ctx.coef_index.unwrap_or(n_coef.saturating_sub(1));
413 replace_coef_moments(&mut prepared.prior, idx, mean, var)?;
414 extra.push_str(&format!(
415 "; reweighted mean={mean:.6} var={var:.6} ess={:.3}",
416 adj.kish_ess()
417 ));
418 }
419 (
420 TransportPolicy::InvariantConditionalOutcome
421 | TransportPolicy::InvariantEffectModifiers,
422 None,
423 ) => {
424 }
426 }
427
428 prepared.prior.restrictions.push(transport_assumption(
429 policy,
430 src.id.as_ref(),
431 sp,
432 ctx.target_population,
433 &extra,
434 ));
435
436 if let Some(a) = alpha_override {
437 prepared.weight.alpha = a;
438 }
439 }
440
441 outcomes.push(TransportOutcome {
442 source_id: Arc::clone(&src.id),
443 required,
444 alpha_override,
445 zero_reason,
446 });
447 out_sources.push(prepared);
448 }
449
450 Ok((out_sources, outcomes))
451}
452
453pub fn compose_with_transport(
462 sources: &[ExternalPriorSource],
463 baseline: &PriorSet,
464 ctx: &TransportContext<'_>,
465) -> Result<(ComposedPrior, Vec<TransportOutcome>), ProbError> {
466 let (prepared, outcomes) = apply_transport(sources, ctx)?;
467 let requested: Vec<f64> = sources.iter().map(|s| s.weight.alpha).collect();
468 let applied: Vec<f64> = prepared
469 .iter()
470 .zip(outcomes.iter())
471 .map(|(s, o)| o.alpha_override.unwrap_or(s.weight.alpha))
472 .collect();
473 let prepared: Vec<ExternalPriorSource> = prepared
475 .into_iter()
476 .zip(applied.iter())
477 .map(|(mut s, &a)| {
478 s.weight.alpha = a;
479 s
480 })
481 .collect();
482 let composed = compose_external_priors_with_alphas(&prepared, &requested, &applied, baseline)?;
483 Ok((composed, outcomes))
484}
485
486#[cfg(test)]
487mod tests {
488 use super::*;
489 use crate::external_prior::ExternalPriorWeight;
490
491 fn gauss(mean: f64, var: f64) -> PriorSet {
492 let mut p = PriorSet::new();
493 p.push(PriorSpec::GaussianCoefficients(
494 GaussianCoefficientPrior::shared(1, mean, var).unwrap(),
495 ));
496 p
497 }
498
499 fn source(id: &str, mean: f64, alpha: f64) -> ExternalPriorSource {
500 ExternalPriorSource {
501 id: Arc::from(id),
502 prior: gauss(mean, 1.0),
503 weight: ExternalPriorWeight::power(alpha).unwrap(),
504 ess: None,
505 }
506 }
507
508 #[test]
509 fn same_population_skips_transport() {
510 let sources = [source("a", 1.0, 0.8)];
511 let ctx = TransportContext {
512 source_populations: &[Some("us")],
513 target_population: Some("us"),
514 policy: None,
515 adjustment: None,
516 coef_index: None,
517 };
518 let (out, outcomes) = apply_transport(&sources, &ctx).unwrap();
519 assert!(!outcomes[0].required);
520 assert!(out[0].prior.restrictions.is_empty());
521 assert!((out[0].weight.alpha - 0.8).abs() < 1e-12);
522 }
523
524 #[test]
525 fn mismatch_without_policy_errors() {
526 let sources = [source("a", 1.0, 1.0)];
527 let ctx = TransportContext {
528 source_populations: &[Some("us")],
529 target_population: Some("eu"),
530 policy: None,
531 adjustment: None,
532 coef_index: None,
533 };
534 let err = apply_transport(&sources, &ctx).unwrap_err();
535 assert_eq!(err.code(), "transport_policy_required");
536 }
537
538 #[test]
539 fn claim_only_records_assumption() {
540 let sources = [source("a", 2.0, 1.0)];
541 let baseline = gauss(0.0, 4.0);
542 let ctx = TransportContext {
543 source_populations: &[Some("us")],
544 target_population: Some("eu"),
545 policy: Some(TransportPolicy::InvariantConditionalOutcome),
546 adjustment: None,
547 coef_index: None,
548 };
549 let (composed, outcomes) = compose_with_transport(&sources, &baseline, &ctx).unwrap();
550 assert!(outcomes[0].required);
551 assert!(outcomes[0].alpha_override.is_none());
552 assert!(
553 composed.prior.restrictions.iter().any(|r| r.id.as_ref() == TRANSPORT_ASSUMPTION_ID)
554 );
555 assert!((composed.alphas_applied[0] - 1.0).abs() < 1e-12);
556 let coef = composed.prior.gaussian_coefficients().unwrap();
557 assert!(coef.mean[0].is_finite());
558 assert!(coef.variance[0].is_finite() && coef.variance[0] > 0.0);
559 }
560
561 #[test]
562 fn propensity_without_weights_forces_alpha_zero() {
563 let sources = [source("a", 2.0, 0.9)];
564 let baseline = gauss(0.0, 4.0);
565 let ctx = TransportContext {
566 source_populations: &[Some("us")],
567 target_population: Some("eu"),
568 policy: Some(TransportPolicy::InvariantPropensity),
569 adjustment: None,
570 coef_index: None,
571 };
572 let (composed, outcomes) = compose_with_transport(&sources, &baseline, &ctx).unwrap();
573 assert_eq!(outcomes[0].alpha_override, Some(0.0));
574 assert!((composed.alphas_requested[0] - 0.9).abs() < 1e-12);
575 assert!((composed.alphas_applied[0] - 0.0).abs() < 1e-12);
576 assert!(
577 composed.prior.restrictions.iter().any(|r| r.id.as_ref() == TRANSPORT_ASSUMPTION_ID)
578 );
579 }
580
581 #[test]
582 fn weighted_moments_shift_mean() {
583 let adj = TransportAdjustment::new([0.0, 0.0, 10.0], [0.0, 0.0, 1.0]).unwrap();
585 let (mean, var) = adj.weighted_moments();
586 assert!((mean - 10.0).abs() < 1e-12);
587 assert!(var >= REWEIGHT_VAR_FLOOR);
588
589 let sources = [source("a", 0.0, 1.0)];
590 let baseline = gauss(0.0, 100.0);
591 let ctx = TransportContext {
592 source_populations: &[Some("us")],
593 target_population: Some("eu"),
594 policy: Some(TransportPolicy::InvariantConditionalOutcome),
595 adjustment: Some(&adj),
596 coef_index: Some(0),
597 };
598 let (composed, _) = compose_with_transport(&sources, &baseline, &ctx).unwrap();
599 let coef = composed.prior.gaussian_coefficients().unwrap();
600 assert!(coef.mean[0] > 5.0, "mean {}", coef.mean[0]);
602 }
603
604 #[test]
605 fn rejects_invalid_weights() {
606 assert!(TransportAdjustment::new([1.0], [-0.1]).is_err());
607 assert!(TransportAdjustment::new([1.0, 2.0], [1.0]).is_err());
608 assert!(TransportAdjustment::new([1.0], [0.0]).is_err());
609 }
610
611 #[test]
612 fn untagged_both_sides_ok_without_policy() {
613 let sources = [source("a", 1.0, 1.0)];
614 let ctx = TransportContext {
615 source_populations: &[None],
616 target_population: None,
617 policy: None,
618 adjustment: None,
619 coef_index: None,
620 };
621 assert!(apply_transport(&sources, &ctx).is_ok());
622 }
623}