Skip to main content

phasesmith_workflows/
quantitative.rs

1//! Quantitative crystalline phase fractions from compatible Rietveld scales.
2
3use std::error::Error;
4use std::fmt::{Display, Formatter};
5
6/// Metadata required by the Hill--Howard scale relation.
7#[derive(Clone, Debug, PartialEq)]
8pub struct QuantitativePhase {
9    /// Stable phase identity.
10    pub phase_id: String,
11    /// Compatible non-negative refined phase scale.
12    pub scale: f64,
13    /// Formula units per crystallographic unit cell.
14    pub formula_units_per_cell: f64,
15    /// Formula mass in grams per mole.
16    pub formula_mass_g_mol: f64,
17    /// Unit-cell volume in cubic ångströms.
18    pub cell_volume_angstrom3: f64,
19}
20
21impl QuantitativePhase {
22    /// Validate one quantitative phase record.
23    ///
24    /// # Errors
25    ///
26    /// Returns [`QuantitativeError`] for empty identity, negative scale, or
27    /// non-positive/non-finite physical metadata.
28    pub fn new(
29        phase_id: impl Into<String>,
30        scale: f64,
31        formula_units_per_cell: f64,
32        formula_mass_g_mol: f64,
33        cell_volume_angstrom3: f64,
34    ) -> Result<Self, QuantitativeError> {
35        let result = Self {
36            phase_id: phase_id.into(),
37            scale,
38            formula_units_per_cell,
39            formula_mass_g_mol,
40            cell_volume_angstrom3,
41        };
42        result.validate()?;
43        Ok(result)
44    }
45
46    fn validate(&self) -> Result<(), QuantitativeError> {
47        if self.phase_id.is_empty() {
48            return Err(QuantitativeError::EmptyPhaseId);
49        }
50        if !self.scale.is_finite() || self.scale < 0.0 {
51            return Err(QuantitativeError::InvalidScale);
52        }
53        if [
54            self.formula_units_per_cell,
55            self.formula_mass_g_mol,
56            self.cell_volume_angstrom3,
57        ]
58        .into_iter()
59        .any(|value| !value.is_finite() || value <= 0.0)
60        {
61            return Err(QuantitativeError::InvalidMetadata);
62        }
63        Ok(())
64    }
65}
66
67/// One normalized crystalline weight fraction.
68#[derive(Clone, Debug, PartialEq)]
69pub struct PhaseWeightFraction {
70    /// Stable phase identity.
71    pub phase_id: String,
72    /// Fraction of the supplied crystalline phases in `[0, 1]`.
73    pub weight_fraction: f64,
74}
75
76/// Weight fractions and their row-major covariance propagated from phase scales.
77#[derive(Clone, Debug, PartialEq)]
78pub struct QuantitativePhaseAnalysis {
79    /// Normalized crystalline weight fractions in input order.
80    pub phases: Vec<PhaseWeightFraction>,
81    /// Row-major fraction covariance with dimension `phases.len()`.
82    pub covariance: Vec<f64>,
83}
84
85/// Calculate Hill--Howard crystalline weight fractions in input order.
86///
87/// `W_p = S_p (Z M V)_p / sum_i S_i (Z M V)_i`.
88///
89/// # Errors
90///
91/// Returns [`QuantitativeError`] for empty, duplicate, invalid, overflowing,
92/// or all-zero phase records.
93pub fn quantitative_phase_analysis(
94    phases: &[QuantitativePhase],
95) -> Result<Vec<PhaseWeightFraction>, QuantitativeError> {
96    if phases.is_empty() {
97        return Err(QuantitativeError::EmptyPhases);
98    }
99    let mut identities = std::collections::BTreeSet::new();
100    let mut contributions = Vec::with_capacity(phases.len());
101    for phase in phases {
102        phase.validate()?;
103        if !identities.insert(&phase.phase_id) {
104            return Err(QuantitativeError::DuplicatePhaseId);
105        }
106        let contribution = phase.scale
107            * phase.formula_units_per_cell
108            * phase.formula_mass_g_mol
109            * phase.cell_volume_angstrom3;
110        if !contribution.is_finite() {
111            return Err(QuantitativeError::ContributionOverflow);
112        }
113        contributions.push(contribution);
114    }
115    let total = contributions.iter().sum::<f64>();
116    if !total.is_finite() || total <= 0.0 {
117        return Err(QuantitativeError::ZeroTotal);
118    }
119    Ok(phases
120        .iter()
121        .zip(contributions)
122        .map(|(phase, contribution)| PhaseWeightFraction {
123            phase_id: phase.phase_id.clone(),
124            weight_fraction: contribution / total,
125        })
126        .collect())
127}
128
129/// Calculate weight fractions and analytically propagate a phase-scale covariance.
130///
131/// For `c_i = S_i (Z M V)_i`, `W_i = c_i / sum(c)`, the scale derivative is
132/// `dW_i/dS_j = (delta_ij k_i - W_i k_j) / sum(c)`, where `k_i = (Z M V)_i`.
133///
134/// # Errors
135///
136/// Returns [`QuantitativeError`] for invalid phases or covariance shape/values.
137pub fn quantitative_phase_analysis_with_covariance(
138    phases: &[QuantitativePhase],
139    scale_covariance: &[f64],
140) -> Result<QuantitativePhaseAnalysis, QuantitativeError> {
141    let fractions = quantitative_phase_analysis(phases)?;
142    let count = phases.len();
143    if scale_covariance.len()
144        != count
145            .checked_mul(count)
146            .ok_or(QuantitativeError::CovarianceShape)?
147    {
148        return Err(QuantitativeError::CovarianceShape);
149    }
150    if scale_covariance.iter().any(|value| !value.is_finite()) {
151        return Err(QuantitativeError::NonFiniteCovariance);
152    }
153    let factors = phases
154        .iter()
155        .map(|phase| {
156            phase.formula_units_per_cell * phase.formula_mass_g_mol * phase.cell_volume_angstrom3
157        })
158        .collect::<Vec<_>>();
159    let total = phases
160        .iter()
161        .zip(&factors)
162        .map(|(phase, factor)| phase.scale * factor)
163        .sum::<f64>();
164    let mut jacobian = vec![0.0; count * count];
165    for row in 0..count {
166        for column in 0..count {
167            jacobian[row * count + column] = (if row == column { factors[row] } else { 0.0 }
168                - fractions[row].weight_fraction * factors[column])
169                / total;
170        }
171    }
172    let mut covariance = vec![0.0; count * count];
173    for row in 0..count {
174        for column in 0..count {
175            let mut value = 0.0;
176            for left in 0..count {
177                for right in 0..count {
178                    value += jacobian[row * count + left]
179                        * scale_covariance[left * count + right]
180                        * jacobian[column * count + right];
181                }
182            }
183            if !value.is_finite() {
184                return Err(QuantitativeError::NonFiniteCovariance);
185            }
186            covariance[row * count + column] = value;
187        }
188    }
189    Ok(QuantitativePhaseAnalysis {
190        phases: fractions,
191        covariance,
192    })
193}
194
195/// Invalid quantitative-phase input.
196#[derive(Clone, Copy, Debug, PartialEq, Eq)]
197pub enum QuantitativeError {
198    /// At least one phase is required.
199    EmptyPhases,
200    /// Phase identity must not be empty.
201    EmptyPhaseId,
202    /// Phase identities must be unique.
203    DuplicatePhaseId,
204    /// Scale must be finite and non-negative.
205    InvalidScale,
206    /// Z, mass, and cell volume must be positive and finite.
207    InvalidMetadata,
208    /// A scale-times-metadata product overflowed.
209    ContributionOverflow,
210    /// At least one phase scale must be positive.
211    ZeroTotal,
212    /// Phase-scale covariance did not have square phase dimension.
213    CovarianceShape,
214    /// Phase-scale covariance or propagated covariance was not finite.
215    NonFiniteCovariance,
216}
217
218impl Display for QuantitativeError {
219    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
220        formatter.write_str(match self {
221            Self::EmptyPhases => "quantitative phase analysis requires at least one phase",
222            Self::EmptyPhaseId => "quantitative phase IDs must not be empty",
223            Self::DuplicatePhaseId => "quantitative phase IDs must be unique",
224            Self::InvalidScale => "quantitative phase scales must be finite and non-negative",
225            Self::InvalidMetadata => "quantitative Z, mass, and volume must be positive and finite",
226            Self::ContributionOverflow => "quantitative phase contribution overflowed",
227            Self::ZeroTotal => "at least one quantitative phase scale must be positive",
228            Self::CovarianceShape => "scale covariance must be square with one row per phase",
229            Self::NonFiniteCovariance => {
230                "scale covariance and propagated covariance must be finite"
231            }
232        })
233    }
234}
235
236impl Error for QuantitativeError {}