phasesmith_workflows/
quantitative.rs1use std::error::Error;
4use std::fmt::{Display, Formatter};
5
6#[derive(Clone, Debug, PartialEq)]
8pub struct QuantitativePhase {
9 pub phase_id: String,
11 pub scale: f64,
13 pub formula_units_per_cell: f64,
15 pub formula_mass_g_mol: f64,
17 pub cell_volume_angstrom3: f64,
19}
20
21impl QuantitativePhase {
22 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#[derive(Clone, Debug, PartialEq)]
69pub struct PhaseWeightFraction {
70 pub phase_id: String,
72 pub weight_fraction: f64,
74}
75
76#[derive(Clone, Debug, PartialEq)]
78pub struct QuantitativePhaseAnalysis {
79 pub phases: Vec<PhaseWeightFraction>,
81 pub covariance: Vec<f64>,
83}
84
85pub 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
129pub 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#[derive(Clone, Copy, Debug, PartialEq, Eq)]
197pub enum QuantitativeError {
198 EmptyPhases,
200 EmptyPhaseId,
202 DuplicatePhaseId,
204 InvalidScale,
206 InvalidMetadata,
208 ContributionOverflow,
210 ZeroTotal,
212 CovarianceShape,
214 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 {}