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