openbnct_core/
exposure.rs

1// SPDX-License-Identifier: Apache-2.0
2
3//! Multi-exposure physical-dose accumulation.
4//!
5//! A research scenario may deliver several weighted exposures — fields,
6//! fractions, or repeated runs — onto one voxel grid. The
7//! `openbnct.exposure-plan/0.1.0` contract declares each exposure's bound
8//! dose bundle, multiplicative delivery weight, basis, and optional
9//! duration and boron assumption. Accumulation sums `weight * dose` and
10//! propagates 1-sigma uncertainties under the declared covariance
11//! assumption; the only supported assumption is statistical independence
12//! between exposures, so sigmas add in quadrature. Within one exposure the
13//! component covariance is already accounted for by that bundle's
14//! dedicated physical-total estimator — the accumulated total sums the
15//! exposure totals, not the accumulated components.
16
17use std::collections::BTreeSet;
18
19use serde::{Deserialize, Serialize};
20use thiserror::Error;
21
22use crate::{
23    ContentReference, DoseVolume, PHYSICAL_DOSE_BUNDLE_SCHEMA, PhysicalDoseBundle,
24    PhysicalTotalDoseVolume, TotalUncertaintyMethod, ValidationError,
25};
26
27/// Schema identifier carried by every `ExposurePlan`.
28pub const EXPOSURE_PLAN_SCHEMA: &str = "openbnct.exposure-plan/0.1.0";
29
30/// A file on disk bound by content hash. `path` is resolved relative to the
31/// plan file's directory by the caller; the referenced bytes must hash to
32/// `sha256`.
33#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
34#[serde(deny_unknown_fields)]
35pub struct BoundFileReference {
36    pub id: String,
37    pub sha256: String,
38    pub path: String,
39}
40
41/// What an exposure's delivery weight physically represents.
42#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
43#[serde(rename_all = "snake_case")]
44pub enum WeightBasis {
45    /// Fraction of a prescribed delivery (for example one fraction of a
46    /// multi-fraction course).
47    DeliveredFraction,
48    /// Ratio of delivered source strength or monitor units to the bundle's
49    /// simulated normalization.
50    SourceStrengthScaling,
51    /// Ratio of delivered particle histories to the simulated histories.
52    DeliveredHistories,
53    /// Weight derived outside the above bases; the plan's notes must say how.
54    Manual,
55}
56
57/// Covariance model between exposures. Only independence is supported:
58/// statistically coupled exposures (for example two fractions sharing one
59/// transport run) must be merged into a single exposure first.
60#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
61#[serde(rename_all = "snake_case")]
62pub enum ExposureCovariance {
63    IndependentExposures,
64}
65
66/// One weighted exposure: a bound dose bundle plus its delivery semantics.
67#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
68#[serde(deny_unknown_fields)]
69pub struct Exposure {
70    pub name: String,
71    pub dose_bundle: BoundFileReference,
72    /// Multiplicative delivery scale applied to the bundle's per-source
73    /// dose. Must be finite and non-negative.
74    pub weight: f64,
75    pub weight_basis: WeightBasis,
76    /// Irradiation duration in seconds, when the plan tracks dose rate.
77    pub duration_s: Option<f64>,
78    /// Free-text record of the boron concentration or compound assumed for
79    /// this exposure — for example the ppm and biodistribution the bound
80    /// run's material encoded.
81    pub boron_assumption: Option<String>,
82}
83
84/// A set of weighted exposures accumulated onto one voxel grid.
85#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
86#[serde(deny_unknown_fields)]
87pub struct ExposurePlan {
88    #[serde(deserialize_with = "crate::deserialize_contract_id")]
89    pub schema_version: String,
90    pub id: String,
91    /// Case identity carried by the accumulated bundle. Exposures may come
92    /// from runs with different case ids (for example different field
93    /// directions); the plan's case id is the accumulated scenario's.
94    pub case_id: String,
95    pub covariance: ExposureCovariance,
96    pub exposures: Vec<Exposure>,
97}
98
99impl ExposurePlan {
100    /// Collect every detectable problem with the plan rather than stopping
101    /// at the first, for user-facing diagnostics on malformed plans.
102    /// Semantic errors that need cross-field context (duplicate exposure
103    /// names) are reported alongside field-level errors.
104    pub fn validate_diagnostics(&self) -> Vec<ExposurePlanError> {
105        let mut issues = Vec::new();
106        if !crate::schema_matches(&self.schema_version, EXPOSURE_PLAN_SCHEMA) {
107            issues.push(ExposurePlanError::UnsupportedSchema(
108                self.schema_version.clone(),
109            ));
110            return issues;
111        }
112        for (label, value) in [("id", self.id.as_str()), ("case_id", self.case_id.as_str())] {
113            if value.trim().is_empty() {
114                issues.push(ExposurePlanError::EmptyIdentifier(label));
115            }
116        }
117        if self.exposures.is_empty() {
118            issues.push(ExposurePlanError::NoExposures);
119            return issues;
120        }
121        let mut names = BTreeSet::new();
122        for exposure in &self.exposures {
123            if exposure.name.trim().is_empty() {
124                issues.push(ExposurePlanError::EmptyIdentifier("exposure.name"));
125            } else if !names.insert(exposure.name.as_str()) {
126                issues.push(ExposurePlanError::DuplicateExposure(exposure.name.clone()));
127            }
128            if !exposure.weight.is_finite() || exposure.weight < 0.0 {
129                issues.push(ExposurePlanError::InvalidWeight(exposure.name.clone()));
130            }
131            if let Some(duration) = exposure.duration_s
132                && (!duration.is_finite() || duration <= 0.0)
133            {
134                issues.push(ExposurePlanError::InvalidDuration(exposure.name.clone()));
135            }
136            let reference = &exposure.dose_bundle;
137            if reference.id.trim().is_empty() || reference.path.trim().is_empty() {
138                issues.push(ExposurePlanError::EmptyIdentifier("exposure.dose_bundle"));
139            }
140            if (ContentReference {
141                id: reference.id.clone(),
142                sha256: reference.sha256.clone(),
143            })
144            .validate()
145            .is_err()
146            {
147                issues.push(ExposurePlanError::InvalidBundleReference(
148                    exposure.name.clone(),
149                ));
150            }
151        }
152        issues
153    }
154
155    pub fn validate(&self) -> Result<(), ExposurePlanError> {
156        if !crate::schema_matches(&self.schema_version, EXPOSURE_PLAN_SCHEMA) {
157            return Err(ExposurePlanError::UnsupportedSchema(
158                self.schema_version.clone(),
159            ));
160        }
161        for (label, value) in [("id", self.id.as_str()), ("case_id", self.case_id.as_str())] {
162            if value.trim().is_empty() {
163                return Err(ExposurePlanError::EmptyIdentifier(label));
164            }
165        }
166        if self.exposures.is_empty() {
167            return Err(ExposurePlanError::NoExposures);
168        }
169        let mut names = BTreeSet::new();
170        for exposure in &self.exposures {
171            if exposure.name.trim().is_empty() {
172                return Err(ExposurePlanError::EmptyIdentifier("exposure.name"));
173            }
174            if !names.insert(exposure.name.as_str()) {
175                return Err(ExposurePlanError::DuplicateExposure(exposure.name.clone()));
176            }
177            if !exposure.weight.is_finite() || exposure.weight < 0.0 {
178                return Err(ExposurePlanError::InvalidWeight(exposure.name.clone()));
179            }
180            if let Some(duration) = exposure.duration_s
181                && (!duration.is_finite() || duration <= 0.0)
182            {
183                return Err(ExposurePlanError::InvalidDuration(exposure.name.clone()));
184            }
185            let reference = &exposure.dose_bundle;
186            if reference.id.trim().is_empty() || reference.path.trim().is_empty() {
187                return Err(ExposurePlanError::EmptyIdentifier("exposure.dose_bundle"));
188            }
189            ContentReference {
190                id: reference.id.clone(),
191                sha256: reference.sha256.clone(),
192            }
193            .validate()
194            .map_err(|_| ExposurePlanError::InvalidBundleReference(exposure.name.clone()))?;
195        }
196        Ok(())
197    }
198}
199
200/// Accumulate the plan's exposure bundles into one physical dose bundle.
201///
202/// `bundles` must align one-to-one with `plan.exposures` — the caller
203/// resolves each bound path, verifies the recorded SHA-256, and parses the
204/// bundle before calling. Every bundle must share the accumulated grid's
205/// geometry, component profile, component set, and dose unit. The result
206/// reuses `openbnct.physical-dose-bundle/0.2.0`; when the exposures bind
207/// different response sets (for example differing boron loading), the
208/// output's `response_set` reference points at the plan artifact, which
209/// enumerates the actual bound sets.
210pub fn accumulate_exposures(
211    plan: &ExposurePlan,
212    plan_sha256: &str,
213    bundles: &[PhysicalDoseBundle],
214) -> Result<PhysicalDoseBundle, ExposurePlanError> {
215    plan.validate()?;
216    if bundles.len() != plan.exposures.len() {
217        return Err(ExposurePlanError::BundleCountMismatch {
218            expected: plan.exposures.len(),
219            actual: bundles.len(),
220        });
221    }
222    for bundle in bundles {
223        bundle.validate()?;
224    }
225    let first = &bundles[0];
226    for (exposure, bundle) in plan.exposures.iter().zip(bundles.iter()) {
227        if bundle.geometry != first.geometry {
228            return Err(ExposurePlanError::GeometryMismatch(exposure.name.clone()));
229        }
230        if bundle.component_profile != first.component_profile {
231            return Err(ExposurePlanError::ComponentProfileMismatch(
232                exposure.name.clone(),
233            ));
234        }
235        if bundle.physical_total.unit != first.physical_total.unit {
236            return Err(ExposurePlanError::DoseUnitMismatch(exposure.name.clone()));
237        }
238        if bundle.frame_of_reference_uid != first.frame_of_reference_uid {
239            return Err(ExposurePlanError::GeometryMismatch(exposure.name.clone()));
240        }
241        let first_components: BTreeSet<_> = first.components.iter().map(|v| v.component).collect();
242        let components: BTreeSet<_> = bundle.components.iter().map(|v| v.component).collect();
243        if components != first_components {
244            return Err(ExposurePlanError::ComponentSetMismatch(
245                exposure.name.clone(),
246            ));
247        }
248    }
249    let voxel_count = first.physical_total.values.len();
250    let shared_response_set = if bundles
251        .iter()
252        .all(|bundle| bundle.response_set == first.response_set)
253    {
254        first.response_set.clone()
255    } else {
256        ContentReference {
257            id: format!("multiple:{}", plan.id),
258            sha256: plan_sha256.to_owned(),
259        }
260    };
261
262    let mut components = Vec::with_capacity(first.components.len());
263    for component in &first.components {
264        let mut values = vec![0.0; voxel_count];
265        let mut sigmas = vec![0.0_f64; voxel_count];
266        let mut have_sigmas = true;
267        for (exposure, bundle) in plan.exposures.iter().zip(bundles.iter()) {
268            let volume = bundle
269                .components
270                .iter()
271                .find(|v| v.component == component.component)
272                .expect("component set checked above");
273            for (index, value) in values.iter_mut().enumerate() {
274                *value += exposure.weight * volume.values[index];
275            }
276            match &volume.absolute_standard_uncertainty {
277                Some(sigma) => {
278                    for (index, accumulated) in sigmas.iter_mut().enumerate() {
279                        *accumulated += (exposure.weight * sigma[index]).powi(2);
280                    }
281                }
282                None => have_sigmas = false,
283            }
284        }
285        components.push(DoseVolume {
286            component: component.component,
287            unit: component.unit,
288            values,
289            absolute_standard_uncertainty: have_sigmas
290                .then(|| sigmas.iter().map(|sigma| sigma.sqrt()).collect()),
291        });
292    }
293
294    let mut total_values = vec![0.0; voxel_count];
295    let mut total_sigmas = vec![0.0_f64; voxel_count];
296    let mut total_method = TotalUncertaintyMethod::DedicatedEstimator;
297    let mut have_total_sigmas = true;
298    for (exposure, bundle) in plan.exposures.iter().zip(bundles.iter()) {
299        for (index, value) in total_values.iter_mut().enumerate() {
300            *value += exposure.weight * bundle.physical_total.values[index];
301        }
302        match (
303            &bundle.physical_total.absolute_standard_uncertainty,
304            bundle.physical_total.uncertainty_method,
305        ) {
306            (Some(sigma), TotalUncertaintyMethod::DedicatedEstimator) => {
307                for (index, accumulated) in total_sigmas.iter_mut().enumerate() {
308                    *accumulated += (exposure.weight * sigma[index]).powi(2);
309                }
310            }
311            (Some(sigma), TotalUncertaintyMethod::BatchCovariance) => {
312                total_method = TotalUncertaintyMethod::BatchCovariance;
313                for (index, accumulated) in total_sigmas.iter_mut().enumerate() {
314                    *accumulated += (exposure.weight * sigma[index]).powi(2);
315                }
316            }
317            _ => {
318                total_method = TotalUncertaintyMethod::Unavailable;
319                have_total_sigmas = false;
320            }
321        }
322    }
323    if !have_total_sigmas {
324        total_method = TotalUncertaintyMethod::Unavailable;
325    }
326
327    let bundle = PhysicalDoseBundle {
328        schema_version: PHYSICAL_DOSE_BUNDLE_SCHEMA.into(),
329        case_id: plan.case_id.clone(),
330        frame_of_reference_uid: first.frame_of_reference_uid.clone(),
331        geometry: first.geometry.clone(),
332        component_profile: first.component_profile.clone(),
333        response_set: shared_response_set,
334        components,
335        physical_total: PhysicalTotalDoseVolume {
336            unit: first.physical_total.unit,
337            values: total_values,
338            absolute_standard_uncertainty: have_total_sigmas
339                .then(|| total_sigmas.iter().map(|sigma| sigma.sqrt()).collect()),
340            uncertainty_method: total_method,
341        },
342        provenance_id: format!(
343            "exposure-plan:sha256:{plan_sha256};covariance:independent_exposures"
344        ),
345    };
346    bundle.validate()?;
347    Ok(bundle)
348}
349
350#[derive(Debug, Error)]
351pub enum ExposurePlanError {
352    #[error("unsupported exposure-plan schema {0:?}")]
353    UnsupportedSchema(String),
354    #[error("required identifier {0} is empty")]
355    EmptyIdentifier(&'static str),
356    #[error("exposure plan contains no exposures")]
357    NoExposures,
358    #[error("exposure {0} occurs more than once")]
359    DuplicateExposure(String),
360    #[error("exposure {0} has a non-finite or negative weight")]
361    InvalidWeight(String),
362    #[error("exposure {0} has a non-finite or non-positive duration")]
363    InvalidDuration(String),
364    #[error("exposure {0} carries an invalid dose-bundle content reference")]
365    InvalidBundleReference(String),
366    #[error("plan declares {expected} exposures but {actual} bundles were supplied")]
367    BundleCountMismatch { expected: usize, actual: usize },
368    #[error("exposure {0} bundle geometry or frame of reference differs")]
369    GeometryMismatch(String),
370    #[error("exposure {0} binds a different component profile")]
371    ComponentProfileMismatch(String),
372    #[error("exposure {0} binds a different component set")]
373    ComponentSetMismatch(String),
374    #[error("exposure {0} uses a different dose unit")]
375    DoseUnitMismatch(String),
376    #[error(transparent)]
377    Bundle(#[from] ValidationError),
378}
379
380#[cfg(test)]
381mod tests {
382    use super::*;
383    use crate::{DoseComponent, DoseUnit, GridGeometry};
384
385    fn sha() -> String {
386        "a".repeat(64)
387    }
388
389    fn bundle(scale: f64, sigma_scale: f64) -> PhysicalDoseBundle {
390        let component = |component: DoseComponent, value: f64| DoseVolume {
391            component,
392            unit: DoseUnit::GrayPerSourceParticle,
393            values: vec![value * scale, 2.0 * value * scale],
394            absolute_standard_uncertainty: Some(vec![
395                value * sigma_scale,
396                2.0 * value * sigma_scale,
397            ]),
398        };
399        PhysicalDoseBundle {
400            schema_version: PHYSICAL_DOSE_BUNDLE_SCHEMA.into(),
401            case_id: "case-a".into(),
402            frame_of_reference_uid: Some("1.2.3".into()),
403            geometry: GridGeometry {
404                shape: [2, 1, 1],
405                spacing_mm: [10.0, 10.0, 10.0],
406                origin_mm: [-5.0, 0.0, 0.0],
407                direction: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
408            },
409            component_profile: ContentReference {
410                id: "profile".into(),
411                sha256: sha(),
412            },
413            response_set: ContentReference {
414                id: "response-set".into(),
415                sha256: sha(),
416            },
417            components: vec![
418                component(DoseComponent::Boron, 1.0e-12),
419                component(DoseComponent::Nitrogen, 2.0e-13),
420                component(DoseComponent::Hydrogen, 5.0e-13),
421                component(DoseComponent::Photon, 3.0e-12),
422            ],
423            physical_total: PhysicalTotalDoseVolume {
424                unit: DoseUnit::GrayPerSourceParticle,
425                values: vec![4.8e-12 * scale, 9.6e-12 * scale],
426                absolute_standard_uncertainty: Some(vec![
427                    1.0e-14 * sigma_scale,
428                    2.0e-14 * sigma_scale,
429                ]),
430                uncertainty_method: TotalUncertaintyMethod::DedicatedEstimator,
431            },
432            provenance_id: "test".into(),
433        }
434    }
435
436    fn plan(weights: &[f64]) -> ExposurePlan {
437        ExposurePlan {
438            schema_version: EXPOSURE_PLAN_SCHEMA.into(),
439            id: "openbnct.test.exposure-plan.v1".into(),
440            case_id: "accumulated-case".into(),
441            covariance: ExposureCovariance::IndependentExposures,
442            exposures: weights
443                .iter()
444                .enumerate()
445                .map(|(index, weight)| Exposure {
446                    name: format!("field-{index}"),
447                    dose_bundle: BoundFileReference {
448                        id: format!("bundle-{index}"),
449                        sha256: sha(),
450                        path: format!("bundle-{index}.json"),
451                    },
452                    weight: *weight,
453                    weight_basis: WeightBasis::DeliveredFraction,
454                    duration_s: Some(600.0),
455                    boron_assumption: Some("10 ppm B-10 in CORE".into()),
456                })
457                .collect(),
458        }
459    }
460
461    #[test]
462    fn accumulates_unequal_weights_with_quadrature_uncertainty() {
463        let plan = plan(&[1.0, 0.5]);
464        let bundles = [bundle(1.0, 1.0e-14), bundle(1.0, 1.0e-14)];
465        let accumulated = accumulate_exposures(&plan, &sha(), &bundles).unwrap();
466        assert_eq!(accumulated.case_id, "accumulated-case");
467        let boron = &accumulated.components[0];
468        // 1.0 * 1e-12 + 0.5 * 1e-12 per unit value.
469        assert_eq!(boron.values, vec![1.5e-12, 3.0e-12]);
470        // sqrt(1^2 + 0.25) * sigma -> 1.118x the single-exposure sigma.
471        let sigma = &boron.absolute_standard_uncertainty.as_ref().unwrap()[0];
472        let expected = 1.0e-14 * 1.0e-12 * (1.0_f64 + 0.25).sqrt();
473        assert!((sigma - expected).abs() / expected < 1.0e-12);
474        assert_eq!(
475            accumulated.physical_total.uncertainty_method,
476            TotalUncertaintyMethod::DedicatedEstimator
477        );
478        assert!(accumulated.provenance_id.contains("independent_exposures"));
479    }
480
481    #[test]
482    fn zero_weight_exposure_contributes_nothing() {
483        let plan = plan(&[1.0, 0.0]);
484        let bundles = [bundle(1.0, 1.0e-14), bundle(9.0, 9.0e-14)];
485        let accumulated = accumulate_exposures(&plan, &sha(), &bundles).unwrap();
486        assert_eq!(accumulated.components[0].values, vec![1.0e-12, 2.0e-12]);
487        assert_eq!(accumulated.physical_total.values, vec![4.8e-12, 9.6e-12]);
488    }
489
490    #[test]
491    fn rejects_mismatched_grids_and_components() {
492        let plan = plan(&[1.0, 1.0]);
493        let mut other = bundle(1.0, 1.0e-14);
494        other.geometry.spacing_mm = [5.0, 10.0, 10.0];
495        assert!(matches!(
496            accumulate_exposures(&plan, &sha(), &[bundle(1.0, 1.0e-14), other.clone()]),
497            Err(ExposurePlanError::GeometryMismatch(_))
498        ));
499        let mut fewer = bundle(1.0, 1.0e-14);
500        fewer.components.pop();
501        // A bundle missing a required component fails its own validation
502        // before the cross-exposure component-set check.
503        assert!(matches!(
504            accumulate_exposures(&plan, &sha(), &[bundle(1.0, 1.0e-14), fewer]),
505            Err(ExposurePlanError::Bundle(
506                ValidationError::MissingComponent(DoseComponent::Photon)
507            ))
508        ));
509        let mut wrong_profile = bundle(1.0, 1.0e-14);
510        wrong_profile.component_profile.sha256 = "b".repeat(64);
511        assert!(matches!(
512            accumulate_exposures(&plan, &sha(), &[bundle(1.0, 1.0e-14), wrong_profile]),
513            Err(ExposurePlanError::ComponentProfileMismatch(_))
514        ));
515    }
516
517    #[test]
518    fn plan_validation_rejects_bad_inputs() {
519        assert!(matches!(
520            plan(&[]).validate(),
521            Err(ExposurePlanError::NoExposures)
522        ));
523        assert!(matches!(
524            plan(&[-1.0]).validate(),
525            Err(ExposurePlanError::InvalidWeight(_))
526        ));
527        let mut duplicate = plan(&[1.0, 1.0]);
528        duplicate.exposures[1].name = "field-0".into();
529        assert!(matches!(
530            duplicate.validate(),
531            Err(ExposurePlanError::DuplicateExposure(_))
532        ));
533        let mut bad_schema = plan(&[1.0]);
534        bad_schema.schema_version = "other/0.0.0".into();
535        assert!(matches!(
536            bad_schema.validate(),
537            Err(ExposurePlanError::UnsupportedSchema(_))
538        ));
539    }
540
541    #[test]
542    fn diagnostics_collect_every_issue() {
543        let mut broken = plan(&[1.0, -2.0]);
544        broken.id = "  ".into();
545        broken.exposures[0].name = String::new();
546        broken.exposures[0].dose_bundle.sha256 = "short".into();
547        let issues = broken.validate_diagnostics();
548        // empty id, empty exposure name, bad weight on field-1, bad ref on field-0
549        assert_eq!(issues.len(), 4);
550        assert!(issues.iter().any(|i| matches!(
551            i,
552            ExposurePlanError::EmptyIdentifier(label) if *label == "id"
553        )));
554        assert!(
555            issues
556                .iter()
557                .any(|i| matches!(i, ExposurePlanError::InvalidWeight(name) if name == "field-1"))
558        );
559        assert!(
560            issues
561                .iter()
562                .any(|i| matches!(i, ExposurePlanError::InvalidBundleReference(_)))
563        );
564        // validate() still reports the first problem only.
565        assert!(broken.validate().is_err());
566        assert!(plan(&[1.0]).validate_diagnostics().is_empty());
567    }
568}
openbnct_core/exposure.rs

openbnct_core/
exposure.rs
