1use std::error::Error;
4use std::fmt::{Display, Formatter};
5
6use crate::{
7 ConstraintDerivativeMatrix, DifferentiableBackground, LatticeError, LatticeParameterization,
8 PreparedRietveldObjective, RietveldCalculation, RietveldCalculationOptions, RietveldError,
9 RietveldGeneralParameterError, RietveldInput, RietveldInstrumentParameter,
10 RietveldObjectiveError, RietveldParameterLayout, calculate_rietveld_pattern,
11};
12
13pub const DEFAULT_MAX_LINEARIZATION_ELEMENTS: usize = 10_000_000;
19
20pub struct PreparedGeneralFreeLinearization {
27 calculation: RietveldCalculation,
28 weighted_jacobian: Vec<f64>,
29 sample_scale: Vec<f64>,
30 parameter_count: usize,
31}
32
33impl PreparedGeneralFreeLinearization {
34 #[must_use]
36 pub const fn calculation(&self) -> &RietveldCalculation {
37 &self.calculation
38 }
39
40 #[must_use]
42 pub const fn parameter_count(&self) -> usize {
43 self.parameter_count
44 }
45
46 pub fn jvp(&self, direction: &[f64]) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
53 if direction.len() != self.parameter_count {
54 return Err(RietveldGeneralParameterError::ValueLengthMismatch.into());
55 }
56 Ok(dense_forward_product(
57 &self.weighted_jacobian,
58 direction,
59 self.sample_scale.len(),
60 ))
61 }
62
63 pub fn vjp(&self, samples: &[f64]) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
70 if samples.len() != self.sample_scale.len() {
71 return Err(RietveldGeneralObjectiveError::SampleLengthMismatch);
72 }
73 Ok(dense_reverse_product(
74 &self.weighted_jacobian,
75 samples,
76 self.parameter_count,
77 ))
78 }
79
80 pub fn normal_product(
86 &self,
87 direction: &[f64],
88 damping: f64,
89 ) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
90 if !damping.is_finite() || damping < 0.0 {
91 return Err(RietveldGeneralObjectiveError::InvalidDamping);
92 }
93 let product = self.jvp(direction)?;
94 let mut result = self.vjp(&product)?;
95 for (value, direction) in result.iter_mut().zip(direction) {
96 *value += damping * direction;
97 }
98 Ok(result)
99 }
100
101 pub fn gradient(&self, observed: &[f64]) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
108 if observed.len() != self.sample_scale.len() {
109 return Err(RietveldGeneralObjectiveError::SampleLengthMismatch);
110 }
111 let weighted_residual = self
112 .calculation
113 .y
114 .iter()
115 .zip(observed)
116 .zip(&self.sample_scale)
117 .map(|((calculated, observed), scale)| (calculated - observed) * scale)
118 .collect::<Vec<_>>();
119 self.vjp(&weighted_residual)
120 }
121}
122
123pub struct PreparedGeneralRietveldObjective {
125 input: RietveldInput,
126 options: RietveldCalculationOptions,
127 layout: RietveldParameterLayout,
128 structural: PreparedRietveldObjective,
129 calculation: RietveldCalculation,
130 dense_structural_jacobian: Option<Vec<f64>>,
131 explicit_columns: Vec<(usize, Vec<f64>)>,
132}
133
134impl PreparedGeneralRietveldObjective {
135 pub fn new(
146 input: RietveldInput,
147 options: RietveldCalculationOptions,
148 layout: RietveldParameterLayout,
149 ) -> Result<Self, RietveldGeneralObjectiveError> {
150 Self::new_with_max_linearization_elements(
151 input,
152 options,
153 layout,
154 DEFAULT_MAX_LINEARIZATION_ELEMENTS,
155 )
156 }
157
158 pub fn new_with_max_linearization_elements(
169 input: RietveldInput,
170 options: RietveldCalculationOptions,
171 layout: RietveldParameterLayout,
172 max_linearization_elements: usize,
173 ) -> Result<Self, RietveldGeneralObjectiveError> {
174 let structural = PreparedRietveldObjective::new(
175 input.clone(),
176 options.clone(),
177 layout.structural_layout().clone(),
178 )?;
179 let dense_enabled = max_linearization_elements > 0
180 && structural.dense_element_count()? <= max_linearization_elements;
181 let (calculation, dense_structural_jacobian) = if dense_enabled {
182 let linearization = structural.linearize()?;
183 (linearization.calculation, Some(linearization.jacobian))
184 } else {
185 (calculate_rietveld_pattern(&input, &options)?, None)
186 };
187 let mut explicit_columns = instrument_columns(&input, &calculation, &layout)?;
188 append_sample_physics_columns(&input, &calculation, &layout, &mut explicit_columns)?;
189 append_background_columns(&input, &layout, &mut explicit_columns)?;
190 Ok(Self {
191 input,
192 options,
193 layout,
194 structural,
195 calculation,
196 dense_structural_jacobian,
197 explicit_columns,
198 })
199 }
200
201 #[must_use]
203 pub const fn uses_dense_linearization(&self) -> bool {
204 self.dense_structural_jacobian.is_some()
205 }
206
207 #[must_use]
209 pub const fn preparation_evaluation_count(&self) -> usize {
210 if self.uses_dense_linearization() {
211 1
212 } else {
213 2
214 }
215 }
216
217 #[must_use]
219 pub const fn normal_product_evaluation_count(&self) -> usize {
220 if self.uses_dense_linearization() {
221 0
222 } else {
223 2
224 }
225 }
226
227 #[must_use]
229 pub const fn jvp_evaluation_count(&self) -> usize {
230 if self.uses_dense_linearization() {
231 0
232 } else {
233 1
234 }
235 }
236
237 #[must_use]
239 pub const fn layout(&self) -> &RietveldParameterLayout {
240 &self.layout
241 }
242
243 #[must_use]
245 pub const fn calculation(&self) -> &RietveldCalculation {
246 &self.calculation
247 }
248
249 pub fn free_linearization(
259 &self,
260 derivative: &ConstraintDerivativeMatrix,
261 ) -> Result<Option<PreparedGeneralFreeLinearization>, RietveldGeneralObjectiveError> {
262 if !self.uses_dense_linearization() {
263 return Ok(None);
264 }
265 if derivative.rows != self.layout.parameters().specs().len() {
266 return Err(RietveldGeneralParameterError::ValueLengthMismatch.into());
267 }
268 let sample_count = self.input.pattern.sample_count();
269 let sample_scale = self.sample_scale();
270 let element_count = derivative
271 .columns
272 .checked_mul(sample_count)
273 .ok_or(RietveldObjectiveError::AllocationOverflow)?;
274 let mut weighted_jacobian = Vec::with_capacity(element_count);
275 for free in 0..derivative.columns {
276 let physical = derivative
277 .values
278 .chunks_exact(derivative.columns)
279 .map(|row| row[free])
280 .collect::<Vec<_>>();
281 let (_, mut row) = self.jvp(&physical)?;
282 for (value, scale) in row.iter_mut().zip(&sample_scale) {
283 *value *= scale;
284 }
285 weighted_jacobian.extend(row);
286 }
287 Ok(Some(PreparedGeneralFreeLinearization {
288 calculation: self.calculation.clone(),
289 weighted_jacobian,
290 sample_scale,
291 parameter_count: derivative.columns,
292 }))
293 }
294
295 pub fn jvp(
302 &self,
303 direction: &[f64],
304 ) -> Result<(Vec<f64>, Vec<f64>), RietveldGeneralObjectiveError> {
305 if direction.len() != self.layout.parameters().specs().len() {
306 return Err(RietveldGeneralParameterError::ValueLengthMismatch.into());
307 }
308 let structural_direction = self.layout.structural_direction(direction)?;
309 let (profile, mut derivative) = if let Some(jacobian) = &self.dense_structural_jacobian {
310 (
311 self.calculation.profile_y.clone(),
312 dense_forward_product(
313 jacobian,
314 &structural_direction,
315 self.input.pattern.sample_count(),
316 ),
317 )
318 } else {
319 self.structural.jvp(&structural_direction)?
320 };
321 for (parameter, column) in &self.explicit_columns {
322 let coefficient = direction[*parameter];
323 for (target, value) in derivative.iter_mut().zip(column) {
324 *target += coefficient * value;
325 }
326 }
327 Ok((profile, derivative))
328 }
329
330 pub fn vjp(&self, sample_weights: &[f64]) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
337 if sample_weights.len() != self.input.pattern.sample_count() {
338 return Err(RietveldGeneralObjectiveError::SampleLengthMismatch);
339 }
340 let structural = if let Some(jacobian) = &self.dense_structural_jacobian {
341 dense_reverse_product(
342 jacobian,
343 sample_weights,
344 self.layout.structural_layout().parameters().specs().len(),
345 )
346 } else {
347 self.structural.vjp(sample_weights)?
348 };
349 let mut result = self.layout.expand_structural_gradient(&structural)?;
350 for (parameter, column) in &self.explicit_columns {
351 result[*parameter] += column
352 .iter()
353 .zip(sample_weights)
354 .map(|(left, right)| left * right)
355 .sum::<f64>();
356 }
357 Ok(result)
358 }
359
360 pub fn normal_product(
367 &self,
368 direction: &[f64],
369 damping: f64,
370 ) -> Result<Vec<f64>, RietveldGeneralObjectiveError> {
371 if !damping.is_finite() || damping < 0.0 {
372 return Err(RietveldGeneralObjectiveError::InvalidDamping);
373 }
374 let (_, derivative) = self.jvp(direction)?;
375 let weighted = self.weight_samples(&derivative);
376 let mut result = self.vjp(&weighted)?;
377 for (value, direction) in result.iter_mut().zip(direction) {
378 *value += damping * direction;
379 }
380 Ok(result)
381 }
382
383 pub fn gradient(&self) -> Result<(Vec<f64>, Vec<f64>), RietveldGeneralObjectiveError> {
389 let observed = self
390 .input
391 .pattern
392 .observed_y
393 .as_ref()
394 .ok_or(RietveldError::MissingObservations)?;
395 let residual = self
396 .calculation
397 .y
398 .iter()
399 .zip(observed)
400 .map(|(calculated, observed)| calculated - observed)
401 .collect::<Vec<_>>();
402 Ok((
403 self.calculation.y.clone(),
404 self.vjp(&self.weight_samples(&residual))?,
405 ))
406 }
407
408 fn weight_samples(&self, values: &[f64]) -> Vec<f64> {
409 let mask = self.input.pattern.mask.as_deref();
410 let uncertainty = self
411 .options
412 .use_uncertainty
413 .then_some(self.input.pattern.uncertainty.as_deref())
414 .flatten();
415 values
416 .iter()
417 .enumerate()
418 .map(|(index, value)| {
419 if mask.is_some_and(|mask| !mask[index]) {
420 0.0
421 } else if let Some(sigma) = uncertainty {
422 value / (sigma[index] * sigma[index])
423 } else {
424 *value
425 }
426 })
427 .collect()
428 }
429
430 fn sample_scale(&self) -> Vec<f64> {
431 let mask = self.input.pattern.mask.as_deref();
432 let uncertainty = self
433 .options
434 .use_uncertainty
435 .then_some(self.input.pattern.uncertainty.as_deref())
436 .flatten();
437 (0..self.input.pattern.sample_count())
438 .map(|index| {
439 if mask.is_some_and(|mask| !mask[index]) {
440 0.0
441 } else {
442 uncertainty.map_or(1.0, |sigma| sigma[index].recip())
443 }
444 })
445 .collect()
446 }
447}
448
449fn dense_forward_product(jacobian: &[f64], direction: &[f64], sample_count: usize) -> Vec<f64> {
450 let mut result = vec![0.0; sample_count];
451 for (coefficient, row) in direction.iter().zip(jacobian.chunks_exact(sample_count)) {
452 if *coefficient == 0.0 {
453 continue;
454 }
455 for (target, value) in result.iter_mut().zip(row) {
456 *target += coefficient * value;
457 }
458 }
459 result
460}
461
462fn dense_reverse_product(
463 jacobian: &[f64],
464 sample_weights: &[f64],
465 parameter_count: usize,
466) -> Vec<f64> {
467 jacobian
468 .chunks_exact(sample_weights.len())
469 .take(parameter_count)
470 .map(|row| {
471 row.iter()
472 .zip(sample_weights)
473 .map(|(left, right)| left * right)
474 .sum()
475 })
476 .collect()
477}
478
479fn instrument_columns(
480 input: &RietveldInput,
481 calculation: &RietveldCalculation,
482 layout: &RietveldParameterLayout,
483) -> Result<Vec<(usize, Vec<f64>)>, RietveldGeneralObjectiveError> {
484 let sample_count = input.pattern.sample_count();
485 let mut result = Vec::new();
486 for (parameter, index) in layout.instrument_indices() {
487 let row = instrument_global_row(input, *parameter)?;
488 let mut column = vec![0.0; sample_count];
489 for phase in &calculation.phases {
490 let global = phase
491 .result
492 .accumulation
493 .derivatives
494 .global
495 .as_ref()
496 .ok_or(RietveldGeneralObjectiveError::MissingGlobalDerivatives)?;
497 if row >= global.parameter_count || global.sample_count != sample_count {
498 return Err(RietveldGeneralObjectiveError::GlobalDerivativeShape);
499 }
500 let values = global
501 .values
502 .get(row * sample_count..(row + 1) * sample_count)
503 .ok_or(RietveldGeneralObjectiveError::GlobalDerivativeShape)?;
504 for (target, value) in column.iter_mut().zip(values) {
505 *target += value;
506 }
507 }
508 result.push((*index, column));
509 }
510 Ok(result)
511}
512
513fn append_sample_physics_columns(
514 input: &RietveldInput,
515 calculation: &RietveldCalculation,
516 layout: &RietveldParameterLayout,
517 explicit_columns: &mut Vec<(usize, Vec<f64>)>,
518) -> Result<(), RietveldGeneralObjectiveError> {
519 let mut terms = layout
520 .sample_physics_indices()
521 .map(|(phase, name, parameter)| (phase, name.to_owned(), parameter, 1.0))
522 .collect::<Vec<_>>();
523 for (phase_index, phase) in input.phases.iter().enumerate() {
524 let (_, names) =
525 phase.resolved_sample_physics(input.instrument, input.position_correction)?;
526 if !names
527 .iter()
528 .any(|name| name.starts_with("march_dollase.cell."))
529 {
530 continue;
531 }
532 let parameterization = LatticeParameterization::new(
533 phase.definition().space_group.clone(),
534 phase.definition().cell,
535 )?;
536 let lattice_values = parameterization.values_from_cell(phase.definition().cell)?;
537 let jacobian = parameterization.cell_jacobian(&lattice_values)?;
538 let columns = lattice_values.len();
539 for (parameter_index, spec) in layout.parameters().specs().iter().enumerate() {
540 if spec.key().module() != "lattice"
541 || spec.key().owner_id() != phase.phase_id().as_str()
542 {
543 continue;
544 }
545 let column = parameterization
546 .parameter_names()
547 .iter()
548 .position(|name| name == spec.key().name())
549 .ok_or(RietveldGeneralObjectiveError::GlobalDerivativeShape)?;
550 for (cell_row, name) in [
551 "a_angstrom",
552 "b_angstrom",
553 "c_angstrom",
554 "alpha_deg",
555 "beta_deg",
556 "gamma_deg",
557 ]
558 .iter()
559 .enumerate()
560 {
561 let coefficient = jacobian[cell_row * columns + column];
562 if coefficient != 0.0 {
563 terms.push((
564 phase_index,
565 format!("march_dollase.cell.{name}"),
566 parameter_index,
567 coefficient,
568 ));
569 }
570 }
571 }
572 }
573 for (phase_index, name, parameter_index, coefficient) in terms {
574 append_sample_physics_column(
575 input,
576 calculation,
577 explicit_columns,
578 phase_index,
579 &name,
580 parameter_index,
581 coefficient,
582 )?;
583 }
584 Ok(())
585}
586
587#[allow(clippy::too_many_arguments)]
588fn append_sample_physics_column(
589 input: &RietveldInput,
590 calculation: &RietveldCalculation,
591 explicit_columns: &mut Vec<(usize, Vec<f64>)>,
592 phase_index: usize,
593 name: &str,
594 parameter_index: usize,
595 coefficient: f64,
596) -> Result<(), RietveldGeneralObjectiveError> {
597 let phase = input
598 .phases
599 .get(phase_index)
600 .ok_or(RietveldGeneralObjectiveError::GlobalDerivativeShape)?;
601 let (_, names) = phase.resolved_sample_physics(input.instrument, input.position_correction)?;
602 let provider_row = names
603 .iter()
604 .position(|candidate| candidate == name)
605 .ok_or(RietveldGeneralObjectiveError::GlobalDerivativeShape)?;
606 let sample_count = input.pattern.sample_count();
607 let row = sample_physics_global_start(input) + provider_row;
608 let global = calculation.phases[phase_index]
609 .result
610 .accumulation
611 .derivatives
612 .global
613 .as_ref()
614 .ok_or(RietveldGeneralObjectiveError::MissingGlobalDerivatives)?;
615 let values = global
616 .values
617 .get(row * sample_count..(row + 1) * sample_count)
618 .ok_or(RietveldGeneralObjectiveError::GlobalDerivativeShape)?;
619 if let Some((_, column)) = explicit_columns
620 .iter_mut()
621 .find(|(index, _)| *index == parameter_index)
622 {
623 for (target, value) in column.iter_mut().zip(values) {
624 *target += coefficient * value;
625 }
626 } else {
627 explicit_columns.push((
628 parameter_index,
629 values.iter().map(|value| coefficient * value).collect(),
630 ));
631 }
632 Ok(())
633}
634
635fn append_background_columns(
636 input: &RietveldInput,
637 layout: &RietveldParameterLayout,
638 explicit_columns: &mut Vec<(usize, Vec<f64>)>,
639) -> Result<(), RietveldGeneralObjectiveError> {
640 if layout.background_indices().is_empty() {
641 return Ok(());
642 }
643 let background = input
644 .background
645 .as_ref()
646 .ok_or(RietveldGeneralParameterError::MissingBackground)?;
647 let basis = background.basis(&input.pattern.x_deg)?;
648 if basis.columns != layout.background_indices().len()
649 || basis.rows != input.pattern.sample_count()
650 {
651 return Err(RietveldGeneralObjectiveError::BackgroundDerivativeShape);
652 }
653 for (column_index, parameter_index) in layout.background_indices().iter().enumerate() {
654 explicit_columns.push((
655 *parameter_index,
656 basis
657 .column(column_index)
658 .ok_or(RietveldGeneralObjectiveError::BackgroundDerivativeShape)?,
659 ));
660 }
661 Ok(())
662}
663
664fn sample_physics_global_start(input: &RietveldInput) -> usize {
665 7 - usize::from(input.fixed_spectrum.is_some())
666 + usize::from(input.position_correction.bragg_brentano_mm.is_some())
667 + 2 * usize::from(
668 input
669 .position_correction
670 .debye_scherrer_micrometre
671 .is_some(),
672 )
673 + 2 * usize::from(input.axial_geometry.is_some())
674}
675
676fn instrument_global_row(
677 input: &RietveldInput,
678 parameter: RietveldInstrumentParameter,
679) -> Result<usize, RietveldGeneralObjectiveError> {
680 let monochromatic_row = match parameter {
681 RietveldInstrumentParameter::UDeg2 => 0,
682 RietveldInstrumentParameter::VDeg2 => 1,
683 RietveldInstrumentParameter::WDeg2 => 2,
684 RietveldInstrumentParameter::XDeg => 3,
685 RietveldInstrumentParameter::YDeg => 4,
686 RietveldInstrumentParameter::WavelengthAngstrom => 5,
687 RietveldInstrumentParameter::ZeroShiftDeg => 6,
688 RietveldInstrumentParameter::SampleDisplacementMm => {
689 if input.position_correction.bragg_brentano_mm.is_none() {
690 return Err(RietveldGeneralParameterError::InstrumentGeometryMismatch.into());
691 }
692 7
693 }
694 RietveldInstrumentParameter::DisplaceXMicrometre => {
695 if input
696 .position_correction
697 .debye_scherrer_micrometre
698 .is_none()
699 {
700 return Err(RietveldGeneralParameterError::InstrumentGeometryMismatch.into());
701 }
702 7
703 }
704 RietveldInstrumentParameter::DisplaceYMicrometre => {
705 if input
706 .position_correction
707 .debye_scherrer_micrometre
708 .is_none()
709 {
710 return Err(RietveldGeneralParameterError::InstrumentGeometryMismatch.into());
711 }
712 8
713 }
714 };
715 if input.fixed_spectrum.is_some() {
716 if parameter == RietveldInstrumentParameter::WavelengthAngstrom {
717 return Err(RietveldGeneralParameterError::SpectrumWavelengthRefinement.into());
718 }
719 Ok(monochromatic_row - usize::from(monochromatic_row > 5))
720 } else {
721 Ok(monochromatic_row)
722 }
723}
724
725#[derive(Debug)]
727pub enum RietveldGeneralObjectiveError {
728 Parameter(RietveldGeneralParameterError),
730 Structural(RietveldObjectiveError),
732 Rietveld(RietveldError),
734 Lattice(LatticeError),
736 MissingGlobalDerivatives,
738 GlobalDerivativeShape,
740 BackgroundDerivativeShape,
742 SampleLengthMismatch,
744 InvalidDamping,
746}
747
748impl Display for RietveldGeneralObjectiveError {
749 fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
750 match self {
751 Self::Parameter(error) => Display::fmt(error, formatter),
752 Self::Structural(error) => Display::fmt(error, formatter),
753 Self::Rietveld(error) => Display::fmt(error, formatter),
754 Self::Lattice(error) => Display::fmt(error, formatter),
755 Self::MissingGlobalDerivatives => {
756 formatter.write_str("native engine omitted Rietveld global derivatives")
757 }
758 Self::GlobalDerivativeShape => {
759 formatter.write_str("native Rietveld global derivative shape is invalid")
760 }
761 Self::BackgroundDerivativeShape => {
762 formatter.write_str("native Rietveld background derivative shape is invalid")
763 }
764 Self::SampleLengthMismatch => {
765 formatter.write_str("native Rietveld sample weight length is wrong")
766 }
767 Self::InvalidDamping => {
768 formatter.write_str("native Rietveld damping must be finite and non-negative")
769 }
770 }
771 }
772}
773
774impl Error for RietveldGeneralObjectiveError {
775 fn source(&self) -> Option<&(dyn Error + 'static)> {
776 match self {
777 Self::Parameter(error) => Some(error),
778 Self::Structural(error) => Some(error),
779 Self::Rietveld(error) => Some(error),
780 Self::Lattice(error) => Some(error),
781 Self::MissingGlobalDerivatives
782 | Self::GlobalDerivativeShape
783 | Self::BackgroundDerivativeShape
784 | Self::SampleLengthMismatch
785 | Self::InvalidDamping => None,
786 }
787 }
788}
789
790impl From<RietveldGeneralParameterError> for RietveldGeneralObjectiveError {
791 fn from(value: RietveldGeneralParameterError) -> Self {
792 Self::Parameter(value)
793 }
794}
795impl From<RietveldObjectiveError> for RietveldGeneralObjectiveError {
796 fn from(value: RietveldObjectiveError) -> Self {
797 Self::Structural(value)
798 }
799}
800impl From<RietveldError> for RietveldGeneralObjectiveError {
801 fn from(value: RietveldError) -> Self {
802 Self::Rietveld(value)
803 }
804}
805impl From<LatticeError> for RietveldGeneralObjectiveError {
806 fn from(value: LatticeError) -> Self {
807 Self::Lattice(value)
808 }
809}
810impl From<crate::BackgroundError> for RietveldGeneralObjectiveError {
811 fn from(value: crate::BackgroundError) -> Self {
812 Self::Parameter(RietveldGeneralParameterError::Background(value))
813 }
814}