1use std::collections::BTreeSet;
4use std::error::Error;
5use std::fmt::{Display, Formatter};
6
7use nalgebra::{DMatrix, DVector};
8use phasesmith_crystallography::UnitCell;
9use phasesmith_model::RecordId;
10
11use crate::tof_lebail::{initialize_intensities, normal_tof_stop};
12use crate::tof_multibank::{
13 AcceptedBankState, MultiBankCycleCandidate, MultiBankCycleOutcome, aggregate_metrics,
14 calculate_states, evaluate_bank_metrics, fitted_parameter_count, prepare_multibank_cycle,
15 reflection_intensities,
16};
17use crate::{
18 DiagnosticValue, LatticeBounds, LatticeError, LatticeParameterization, RefinementEventKind,
19 RefinementLimits, RefinementRuntime, ResidualEvaluation, RuntimeError, TerminationReason,
20 TofLeBailError, TofLeBailOptions, TofMultiBankCheckpointBank, TofMultiBankError,
21 TofMultiBankInput, TofMultiBankMetrics, TofMultiBankResultBank, tof_lattice_geometry,
22};
23
24#[derive(Clone, Debug, PartialEq)]
26pub struct TofSharedLatticePhase {
27 phase_id: RecordId,
28 parameterization: LatticeParameterization,
29 bounds: LatticeBounds,
30 initial_cell: UnitCell,
31}
32
33impl TofSharedLatticePhase {
34 pub fn new(
41 phase_id: RecordId,
42 parameterization: LatticeParameterization,
43 bounds: LatticeBounds,
44 initial_cell: UnitCell,
45 ) -> Result<Self, TofMultiBankLatticeError> {
46 let result = Self {
47 phase_id,
48 parameterization,
49 bounds,
50 initial_cell,
51 };
52 result.validate_cell(initial_cell)?;
53 Ok(result)
54 }
55
56 #[must_use]
58 pub const fn phase_id(&self) -> &RecordId {
59 &self.phase_id
60 }
61
62 #[must_use]
64 pub const fn parameterization(&self) -> &LatticeParameterization {
65 &self.parameterization
66 }
67
68 #[must_use]
70 pub const fn bounds(&self) -> &LatticeBounds {
71 &self.bounds
72 }
73
74 #[must_use]
76 pub const fn initial_cell(&self) -> UnitCell {
77 self.initial_cell
78 }
79
80 pub(crate) fn validate_cell(
81 &self,
82 cell: UnitCell,
83 ) -> Result<Vec<f64>, TofMultiBankLatticeError> {
84 if self.bounds.parameter_names() != self.parameterization.parameter_names() {
85 return Err(TofMultiBankLatticeError::InvalidModel(
86 "TOF lattice bounds use a different independent-parameter order",
87 ));
88 }
89 let values = self.parameterization.values_from_cell(cell)?;
90 if values
91 .iter()
92 .zip(self.bounds.lower().iter().zip(self.bounds.upper()))
93 .any(|(value, (lower, upper))| value < lower || value > upper)
94 {
95 return Err(TofMultiBankLatticeError::InvalidModel(
96 "TOF lattice cell lies outside its declared bounds",
97 ));
98 }
99 for corner in self.bounds.corner_values() {
100 self.parameterization.to_cell(&corner)?;
101 }
102 Ok(values)
103 }
104}
105
106#[derive(Clone, Debug, PartialEq)]
108pub struct TofMultiBankLatticeInput {
109 pub multibank: TofMultiBankInput,
111 pub lattice_phases: Vec<TofSharedLatticePhase>,
113}
114
115impl TofMultiBankLatticeInput {
116 pub fn validate(&self) -> Result<(), TofMultiBankLatticeError> {
123 self.multibank.validate()?;
124 if self.lattice_phases.is_empty() {
125 return Err(TofMultiBankLatticeError::InvalidModel(
126 "shared TOF lattice refinement requires at least one phase",
127 ));
128 }
129 let mut phase_ids = BTreeSet::new();
130 for lattice in &self.lattice_phases {
131 if !phase_ids.insert(lattice.phase_id.clone()) {
132 return Err(TofMultiBankLatticeError::InvalidModel(
133 "shared TOF lattice phase IDs must be unique",
134 ));
135 }
136 lattice.validate_cell(lattice.initial_cell)?;
137 for bank in &self.multibank.banks {
138 let phase = bank
139 .input
140 .phases
141 .iter()
142 .find(|phase| phase.phase_id() == lattice.phase_id())
143 .ok_or(TofMultiBankLatticeError::InvalidModel(
144 "shared TOF lattice phase is absent from the bank topology",
145 ))?;
146 let geometry = tof_lattice_geometry(
147 lattice.parameterization(),
148 lattice.initial_cell(),
149 phase.hkl(),
150 bank.input.instrument,
151 )?;
152 if geometry.d_spacing_angstrom != phase.d_spacing_angstrom() {
153 return Err(TofMultiBankLatticeError::InvalidModel(
154 "shared TOF lattice initial cell disagrees with bank d-spacings",
155 ));
156 }
157 }
158 }
159 Ok(())
160 }
161
162 pub(crate) fn parameter_count(&self) -> Result<usize, TofMultiBankLatticeError> {
163 self.lattice_phases
164 .iter()
165 .try_fold(0_usize, |count, phase| {
166 count
167 .checked_add(phase.parameterization.parameter_names().len())
168 .ok_or(TofMultiBankLatticeError::AllocationOverflow)
169 })
170 }
171}
172
173#[derive(Clone, Debug, PartialEq)]
175pub struct TofMultiBankLatticeOptions {
176 pub lebail: TofLeBailOptions,
178 pub lattice_damping: f64,
180 pub max_scaled_lattice_step: f64,
182 pub max_lattice_backtracks: usize,
184}
185
186impl TofMultiBankLatticeOptions {
187 pub fn new(
193 lebail: TofLeBailOptions,
194 lattice_damping: f64,
195 max_scaled_lattice_step: f64,
196 max_lattice_backtracks: usize,
197 ) -> Result<Self, TofMultiBankLatticeError> {
198 let result = Self {
199 lebail,
200 lattice_damping,
201 max_scaled_lattice_step,
202 max_lattice_backtracks,
203 };
204 result.validate()?;
205 Ok(result)
206 }
207
208 pub fn validate(&self) -> Result<(), TofMultiBankLatticeError> {
214 self.lebail.validate()?;
215 if !self.lattice_damping.is_finite()
216 || self.lattice_damping < 0.0
217 || !self.max_scaled_lattice_step.is_finite()
218 || self.max_scaled_lattice_step <= 0.0
219 {
220 return Err(TofMultiBankLatticeError::InvalidOptions);
221 }
222 Ok(())
223 }
224}
225
226#[derive(Clone, Debug, PartialEq)]
228pub struct TofSharedLatticeState {
229 pub phase_id: RecordId,
231 pub cell: UnitCell,
233}
234
235#[derive(Clone, Debug, PartialEq)]
237pub struct TofLatticeParameterChange {
238 pub phase_id: RecordId,
240 pub parameter_name: String,
242 pub before: f64,
244 pub after: f64,
246 pub scaled_change: f64,
248}
249
250#[derive(Clone, Debug, PartialEq)]
252pub struct TofMultiBankLatticeIterationRecord {
253 pub iteration: usize,
255 pub bank_metrics: Vec<ResidualEvaluation>,
257 pub metrics: TofMultiBankMetrics,
259 pub maximum_relative_intensity_change: f64,
261 pub maximum_absolute_background_change: f64,
263 pub scaled_lattice_step_norm: f64,
265 pub lattice_parameter_changes: Vec<TofLatticeParameterChange>,
267}
268
269#[derive(Clone, Debug, PartialEq)]
271pub struct TofMultiBankLatticeCheckpoint {
272 pub completed_iterations: usize,
274 pub banks: Vec<TofMultiBankCheckpointBank>,
276 pub lattice_phases: Vec<TofSharedLatticeState>,
278 pub history: Vec<TofMultiBankLatticeIterationRecord>,
280}
281
282impl TofMultiBankLatticeCheckpoint {
283 pub fn validate_for(
290 &self,
291 input: &TofMultiBankLatticeInput,
292 options: &TofMultiBankLatticeOptions,
293 ) -> Result<(), TofMultiBankLatticeError> {
294 input.validate()?;
295 options.validate()?;
296 if self.completed_iterations != self.history.len()
297 || self.completed_iterations > options.lebail.cycles
298 || self.banks.len() != input.multibank.banks.len()
299 || self.lattice_phases.len() != input.lattice_phases.len()
300 {
301 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
302 "shared TOF lattice checkpoint counts differ from the request",
303 ));
304 }
305 for (saved, model) in self.lattice_phases.iter().zip(&input.lattice_phases) {
306 if &saved.phase_id != model.phase_id() {
307 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
308 "shared TOF lattice checkpoint phase order changed",
309 ));
310 }
311 model.validate_cell(saved.cell)?;
312 }
313 for (saved, original) in self.banks.iter().zip(&input.multibank.banks) {
314 validate_checkpoint_bank(saved, original, input, &self.lattice_phases)?;
315 }
316 if self.history.iter().enumerate().any(|(index, record)| {
317 record.iteration != index + 1
318 || record.bank_metrics.len() != input.multibank.banks.len()
319 }) {
320 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
321 "shared TOF lattice checkpoint history is not contiguous and bank-aligned",
322 ));
323 }
324 Ok(())
325 }
326}
327
328#[derive(Clone, Debug, PartialEq)]
330pub struct TofMultiBankLatticeResult {
331 pub banks: Vec<TofMultiBankResultBank>,
333 pub lattice_phases: Vec<TofSharedLatticeState>,
335 pub metrics: TofMultiBankMetrics,
337 pub history: Vec<TofMultiBankLatticeIterationRecord>,
339 pub termination_reason: TerminationReason,
341 pub checkpoint: TofMultiBankLatticeCheckpoint,
343}
344
345pub fn refine_tof_multibank_lattice(
352 input: &TofMultiBankLatticeInput,
353 options: &TofMultiBankLatticeOptions,
354) -> Result<TofMultiBankLatticeResult, TofMultiBankLatticeError> {
355 input.validate()?;
356 options.validate()?;
357 let evaluations_per_cycle = options
358 .max_lattice_backtracks
359 .checked_add(4)
360 .ok_or(TofMultiBankLatticeError::AllocationOverflow)?;
361 let max_evaluations = options
362 .lebail
363 .cycles
364 .checked_mul(evaluations_per_cycle)
365 .ok_or(TofMultiBankLatticeError::AllocationOverflow)?;
366 let limits = RefinementLimits::new(options.lebail.cycles, max_evaluations, None, 1)?;
367 let mut runtime = RefinementRuntime::new(limits, None)?;
368 refine_tof_multibank_lattice_with_runtime(input, options, None, &mut runtime)
369}
370
371#[allow(clippy::too_many_lines)]
382pub fn refine_tof_multibank_lattice_with_runtime(
383 input: &TofMultiBankLatticeInput,
384 options: &TofMultiBankLatticeOptions,
385 checkpoint: Option<&TofMultiBankLatticeCheckpoint>,
386 runtime: &mut RefinementRuntime<TofMultiBankLatticeCheckpoint>,
387) -> Result<TofMultiBankLatticeResult, TofMultiBankLatticeError> {
388 input.validate()?;
389 options.validate()?;
390 let parameter_count = input.parameter_count()?;
391 let restored = restore_state(input, options, checkpoint)?;
392 let mut states = restored.states;
393 let mut cells = restored.cells;
394 let mut history = restored.history;
395 if let Some(checkpoint) = checkpoint {
396 runtime.resume_accepted(checkpoint.completed_iterations)?;
397 }
398 runtime.emit(
399 RefinementEventKind::Start,
400 "tof_multibank_lattice",
401 "shared-cell multi-bank TOF Le Bail refinement started",
402 vec![
403 (
404 "bank_count".to_owned(),
405 DiagnosticValue::Unsigned(input.multibank.banks.len() as u64),
406 ),
407 (
408 "lattice_parameter_count".to_owned(),
409 DiagnosticValue::Unsigned(parameter_count as u64),
410 ),
411 ],
412 )?;
413 let mut accepted_calculations = None;
414 let mut termination = TerminationReason::MaxIterations;
415 for iteration in restored.first_iteration..=options.lebail.cycles {
416 if let Err(error) = runtime.begin_iteration(iteration) {
417 termination = normal_tof_stop(error)?;
418 break;
419 }
420 let candidate = match prepare_multibank_cycle(
421 &input.multibank,
422 &states,
423 &options.lebail,
424 parameter_count,
425 runtime,
426 )? {
427 MultiBankCycleOutcome::Candidate(candidate) => candidate,
428 MultiBankCycleOutcome::Stopped(reason) => {
429 termination = reason;
430 break;
431 }
432 };
433 let lattice = match lattice_update(input, options, &cells, candidate, runtime)? {
434 LatticeUpdateOutcome::Candidate(candidate) => candidate,
435 LatticeUpdateOutcome::Stopped(reason) => {
436 termination = reason;
437 break;
438 }
439 };
440 states = lattice.cycle.states;
441 cells = lattice.cells;
442 accepted_calculations = Some(lattice.cycle.calculations.clone());
443 history.push(TofMultiBankLatticeIterationRecord {
444 iteration,
445 bank_metrics: lattice.cycle.bank_metrics.clone(),
446 metrics: lattice.cycle.metrics,
447 maximum_relative_intensity_change: lattice.cycle.maximum_relative_intensity_change,
448 maximum_absolute_background_change: lattice.cycle.maximum_absolute_background_change,
449 scaled_lattice_step_norm: lattice.step_norm,
450 lattice_parameter_changes: lattice.changes,
451 });
452 let accepted = checkpoint_from_state(input, &states, &cells, &history);
453 runtime.accept_step(Some(&accepted))?;
454 runtime.emit(
455 RefinementEventKind::Iteration,
456 "tof_multibank_lattice_iteration",
457 "shared-cell multi-bank TOF cycle accepted",
458 vec![
459 (
460 "rwp".to_owned(),
461 DiagnosticValue::Float(lattice.cycle.metrics.rwp),
462 ),
463 (
464 "scaled_lattice_step_norm".to_owned(),
465 DiagnosticValue::Float(lattice.step_norm),
466 ),
467 ],
468 )?;
469 }
470 let calculations = match accepted_calculations {
471 Some(calculations) => calculations,
472 None => calculate_states(&input.multibank, &states, &options.lebail)?,
473 };
474 let bank_metrics =
475 evaluate_bank_metrics(&input.multibank, &states, &calculations, &options.lebail)?;
476 let fitted = fitted_parameter_count(&states)?
477 .checked_add(parameter_count)
478 .ok_or(TofMultiBankLatticeError::AllocationOverflow)?;
479 let metrics = aggregate_metrics(&input.multibank, &calculations, &options.lebail, fitted)?;
480 let checkpoint = checkpoint_from_state(input, &states, &cells, &history);
481 checkpoint.validate_for(input, options)?;
482 let banks = input
483 .multibank
484 .banks
485 .iter()
486 .zip(states)
487 .zip(calculations)
488 .zip(bank_metrics)
489 .map(
490 |(((bank, state), calculation), metrics)| TofMultiBankResultBank {
491 bank_id: bank.bank_id.clone(),
492 calculation,
493 metrics,
494 intensities: reflection_intensities(&state.phases),
495 phases: state.phases,
496 background: state.background,
497 },
498 )
499 .collect();
500 let lattice_phases = lattice_states(input, &cells);
501 runtime.emit(
502 RefinementEventKind::Termination,
503 "tof_multibank_lattice",
504 "shared-cell multi-bank TOF refinement terminated",
505 vec![(
506 "termination_reason".to_owned(),
507 DiagnosticValue::String(termination.as_str().to_owned()),
508 )],
509 )?;
510 Ok(TofMultiBankLatticeResult {
511 banks,
512 lattice_phases,
513 metrics,
514 history,
515 termination_reason: termination,
516 checkpoint,
517 })
518}
519
520struct RestoredState {
521 states: Vec<AcceptedBankState>,
522 cells: Vec<UnitCell>,
523 history: Vec<TofMultiBankLatticeIterationRecord>,
524 first_iteration: usize,
525}
526
527fn restore_state(
528 input: &TofMultiBankLatticeInput,
529 options: &TofMultiBankLatticeOptions,
530 checkpoint: Option<&TofMultiBankLatticeCheckpoint>,
531) -> Result<RestoredState, TofMultiBankLatticeError> {
532 let Some(checkpoint) = checkpoint else {
533 let states = input
534 .multibank
535 .banks
536 .iter()
537 .map(|bank| {
538 Ok(AcceptedBankState {
539 phases: initialize_intensities(&bank.input, &options.lebail)?,
540 background: bank.input.background.clone(),
541 })
542 })
543 .collect::<Result<Vec<_>, TofMultiBankLatticeError>>()?;
544 let cells = input
545 .lattice_phases
546 .iter()
547 .map(TofSharedLatticePhase::initial_cell)
548 .collect::<Vec<_>>();
549 return Ok(RestoredState {
550 states: apply_cells(input, &states, &cells)?,
551 cells,
552 history: Vec::new(),
553 first_iteration: 1,
554 });
555 };
556 checkpoint.validate_for(input, options)?;
557 let cells = checkpoint
558 .lattice_phases
559 .iter()
560 .map(|phase| phase.cell)
561 .collect();
562 Ok(RestoredState {
563 states: checkpoint
564 .banks
565 .iter()
566 .map(|bank| AcceptedBankState {
567 phases: bank.phases.clone(),
568 background: bank.background.clone(),
569 })
570 .collect(),
571 cells,
572 history: checkpoint.history.clone(),
573 first_iteration: checkpoint
574 .completed_iterations
575 .checked_add(1)
576 .ok_or(TofMultiBankLatticeError::AllocationOverflow)?,
577 })
578}
579
580struct LatticeUpdate {
581 cycle: MultiBankCycleCandidate,
582 cells: Vec<UnitCell>,
583 step_norm: f64,
584 changes: Vec<TofLatticeParameterChange>,
585}
586
587enum LatticeUpdateOutcome {
588 Candidate(LatticeUpdate),
589 Stopped(TerminationReason),
590}
591
592#[allow(clippy::too_many_lines)]
593fn lattice_update<C>(
594 input: &TofMultiBankLatticeInput,
595 options: &TofMultiBankLatticeOptions,
596 cells: &[UnitCell],
597 mut cycle: MultiBankCycleCandidate,
598 runtime: &mut RefinementRuntime<C>,
599) -> Result<LatticeUpdateOutcome, TofMultiBankLatticeError> {
600 let packed = pack_lattice(input, cells)?;
601 let (jacobian, residual) =
602 lattice_system(input, cells, &cycle, &packed.scales, &options.lebail)?;
603 let normal = jacobian.transpose() * &jacobian;
604 let rhs = jacobian.transpose() * &residual;
605 let mut regularized = normal;
606 for parameter in 0..regularized.nrows() {
607 regularized[(parameter, parameter)] += options.lattice_damping;
608 }
609 let mut step = if let Some(solution) = regularized.lu().solve(&rhs) {
610 solution
611 } else {
612 jacobian
613 .svd(true, true)
614 .solve(&residual, f64::EPSILON)
615 .map_err(|_| TofMultiBankLatticeError::LinearSolve)?
616 };
617 if step.iter().any(|value| !value.is_finite()) {
618 return Err(TofMultiBankLatticeError::LinearSolve);
619 }
620 for parameter in 0..step.len() {
621 let current = packed.values[parameter] / packed.scales[parameter];
622 let lower = (packed.lower[parameter] / packed.scales[parameter] - current)
623 .max(-options.max_scaled_lattice_step);
624 let upper = (packed.upper[parameter] / packed.scales[parameter] - current)
625 .min(options.max_scaled_lattice_step);
626 step[parameter] = step[parameter].clamp(lower, upper);
627 }
628 if step.norm() == 0.0 {
629 return Ok(LatticeUpdateOutcome::Candidate(LatticeUpdate {
630 cycle,
631 cells: cells.to_vec(),
632 step_norm: 0.0,
633 changes: Vec::new(),
634 }));
635 }
636 let parameter_count = input.parameter_count()?;
637 let local_parameter_count = fitted_parameter_count(&cycle.states)?;
638 let total_parameter_count = local_parameter_count
639 .checked_add(parameter_count)
640 .ok_or(TofMultiBankLatticeError::AllocationOverflow)?;
641 let mut factor = 1.0;
642 for _ in 0..=options.max_lattice_backtracks {
643 let trial_values = packed
644 .values
645 .iter()
646 .zip(&packed.scales)
647 .zip(step.iter())
648 .map(|((value, scale), step)| value + factor * scale * step)
649 .collect::<Vec<_>>();
650 let trial_cells = match unpack_lattice(input, &trial_values) {
651 Ok(cells) => cells,
652 Err(error) if recoverable_trial_error(&error) => {
653 factor *= 0.5;
654 continue;
655 }
656 Err(error) => return Err(error),
657 };
658 let trial_states = match apply_cells(input, &cycle.states, &trial_cells) {
659 Ok(states) => states,
660 Err(error) if recoverable_trial_error(&error) => {
661 factor *= 0.5;
662 continue;
663 }
664 Err(error) => return Err(error),
665 };
666 if let Err(error) = runtime.begin_evaluation() {
667 return Ok(LatticeUpdateOutcome::Stopped(normal_tof_stop(error)?));
668 }
669 let trial_calculations =
670 match calculate_states(&input.multibank, &trial_states, &options.lebail) {
671 Ok(calculations) => calculations,
672 Err(error) if recoverable_multibank_trial_error(&error) => {
673 factor *= 0.5;
674 continue;
675 }
676 Err(error) => return Err(error.into()),
677 };
678 let trial_metrics = aggregate_metrics(
679 &input.multibank,
680 &trial_calculations,
681 &options.lebail,
682 total_parameter_count,
683 )?;
684 if trial_metrics.chi_square < cycle.metrics.chi_square {
685 let trial_bank_metrics = evaluate_bank_metrics(
686 &input.multibank,
687 &trial_states,
688 &trial_calculations,
689 &options.lebail,
690 )?;
691 let changes = lattice_changes(input, &packed.values, &trial_values, &packed.scales)?;
692 cycle.states = trial_states;
693 cycle.calculations = trial_calculations;
694 cycle.bank_metrics = trial_bank_metrics;
695 cycle.metrics = trial_metrics;
696 return Ok(LatticeUpdateOutcome::Candidate(LatticeUpdate {
697 cycle,
698 cells: trial_cells,
699 step_norm: factor * step.norm(),
700 changes,
701 }));
702 }
703 factor *= 0.5;
704 }
705 Ok(LatticeUpdateOutcome::Candidate(LatticeUpdate {
706 cycle,
707 cells: cells.to_vec(),
708 step_norm: 0.0,
709 changes: Vec::new(),
710 }))
711}
712
713fn recoverable_trial_error(error: &TofMultiBankLatticeError) -> bool {
714 matches!(
715 error,
716 TofMultiBankLatticeError::Lattice(_)
717 | TofMultiBankLatticeError::InvalidModel(_)
718 | TofMultiBankLatticeError::MultiBank(TofMultiBankError::LeBail(
719 TofLeBailError::Profile(_) | TofLeBailError::InvalidPhase(_)
720 ))
721 )
722}
723
724fn recoverable_multibank_trial_error(error: &TofMultiBankError) -> bool {
725 matches!(
726 error,
727 TofMultiBankError::LeBail(TofLeBailError::Profile(_) | TofLeBailError::InvalidPhase(_))
728 )
729}
730
731pub(crate) struct PackedLattice {
732 pub(crate) values: Vec<f64>,
733 pub(crate) scales: Vec<f64>,
734 pub(crate) lower: Vec<f64>,
735 pub(crate) upper: Vec<f64>,
736}
737
738pub(crate) fn pack_lattice(
739 input: &TofMultiBankLatticeInput,
740 cells: &[UnitCell],
741) -> Result<PackedLattice, TofMultiBankLatticeError> {
742 if cells.len() != input.lattice_phases.len() {
743 return Err(TofMultiBankLatticeError::InternalInvariant);
744 }
745 let mut values = Vec::new();
746 let mut scales = Vec::new();
747 let mut lower = Vec::new();
748 let mut upper = Vec::new();
749 for (phase, cell) in input.lattice_phases.iter().zip(cells) {
750 let phase_values = phase.validate_cell(*cell)?;
751 scales.extend(phase_values.iter().map(|value| value.abs().max(1.0)));
752 values.extend(phase_values);
753 lower.extend_from_slice(phase.bounds.lower());
754 upper.extend_from_slice(phase.bounds.upper());
755 }
756 Ok(PackedLattice {
757 values,
758 scales,
759 lower,
760 upper,
761 })
762}
763
764pub(crate) fn unpack_lattice(
765 input: &TofMultiBankLatticeInput,
766 values: &[f64],
767) -> Result<Vec<UnitCell>, TofMultiBankLatticeError> {
768 if values.len() != input.parameter_count()? {
769 return Err(TofMultiBankLatticeError::InternalInvariant);
770 }
771 let mut offset = 0;
772 input
773 .lattice_phases
774 .iter()
775 .map(|phase| {
776 let end = offset + phase.parameterization.parameter_names().len();
777 let cell = phase.parameterization.to_cell(&values[offset..end])?;
778 phase.validate_cell(cell)?;
779 offset = end;
780 Ok(cell)
781 })
782 .collect()
783}
784
785pub(crate) fn apply_cells(
786 input: &TofMultiBankLatticeInput,
787 states: &[AcceptedBankState],
788 cells: &[UnitCell],
789) -> Result<Vec<AcceptedBankState>, TofMultiBankLatticeError> {
790 if states.len() != input.multibank.banks.len() || cells.len() != input.lattice_phases.len() {
791 return Err(TofMultiBankLatticeError::InternalInvariant);
792 }
793 input
794 .multibank
795 .banks
796 .iter()
797 .zip(states)
798 .map(|(bank, state)| {
799 let mut phases = state.phases.clone();
800 for (model, cell) in input.lattice_phases.iter().zip(cells) {
801 let index = phases
802 .iter()
803 .position(|phase| phase.phase_id() == model.phase_id())
804 .ok_or(TofMultiBankLatticeError::InternalInvariant)?;
805 let geometry = tof_lattice_geometry(
806 model.parameterization(),
807 *cell,
808 phases[index].hkl(),
809 bank.input.instrument,
810 )?;
811 phases[index] = phases[index].with_d_spacings(&geometry.d_spacing_angstrom)?;
812 }
813 Ok(AcceptedBankState {
814 phases,
815 background: state.background.clone(),
816 })
817 })
818 .collect()
819}
820
821#[allow(clippy::too_many_lines)]
822pub(crate) fn lattice_system(
823 input: &TofMultiBankLatticeInput,
824 cells: &[UnitCell],
825 cycle: &MultiBankCycleCandidate,
826 scales: &[f64],
827 options: &TofLeBailOptions,
828) -> Result<(DMatrix<f64>, DVector<f64>), TofMultiBankLatticeError> {
829 let columns = input.parameter_count()?;
830 if columns != scales.len() {
831 return Err(TofMultiBankLatticeError::InternalInvariant);
832 }
833 let rows = input
834 .multibank
835 .banks
836 .iter()
837 .try_fold(0_usize, |count, bank| {
838 let included = bank.input.pattern.mask.as_ref().map_or_else(
839 || bank.input.pattern.sample_count(),
840 |mask| mask.iter().filter(|value| **value).count(),
841 );
842 count
843 .checked_add(included)
844 .ok_or(TofMultiBankLatticeError::AllocationOverflow)
845 })?;
846 if rows < columns {
847 return Err(TofMultiBankLatticeError::LinearSolve);
848 }
849 let mut selected_jacobian = DMatrix::zeros(rows, columns);
850 let mut selected_residual = DVector::zeros(rows);
851 let mut row_offset = 0;
852 for ((bank, state), calculation) in input
853 .multibank
854 .banks
855 .iter()
856 .zip(&cycle.states)
857 .zip(&cycle.calculations)
858 {
859 let samples = bank.input.pattern.sample_count();
860 let mut physical = DMatrix::<f64>::zeros(samples, columns);
861 let local = &calculation.accumulation.derivatives.local;
862 if local.parameter_count != 2 {
863 return Err(TofMultiBankLatticeError::InternalInvariant);
864 }
865 let mut parameter_offset = 0;
866 for (model, cell) in input.lattice_phases.iter().zip(cells) {
867 let phase_index = state
868 .phases
869 .iter()
870 .position(|phase| phase.phase_id() == model.phase_id())
871 .ok_or(TofMultiBankLatticeError::InternalInvariant)?;
872 let phase = &state.phases[phase_index];
873 let geometry = tof_lattice_geometry(
874 model.parameterization(),
875 *cell,
876 phase.hkl(),
877 bank.input.instrument,
878 )?;
879 let phase_columns = geometry.parameter_names.len();
880 let reflection_offset = calculation.phase_offsets[phase_index];
881 for phase_reflection in 0..phase.reflection_ids().len() {
882 let reflection = reflection_offset + phase_reflection;
883 let begin = local.offsets[reflection];
884 let end = local.offsets[reflection + 1];
885 let start = local.starts[reflection];
886 for parameter in 0..phase_columns {
887 let chain = geometry.d_d_spacing_d_parameters
888 [phase_reflection * phase_columns + parameter]
889 * scales[parameter_offset + parameter];
890 for active in begin..end {
891 physical[(start + active - begin, parameter_offset + parameter)] +=
892 local.values[active * local.parameter_count + 1] * chain;
893 }
894 }
895 }
896 parameter_offset += phase_columns;
897 }
898 if parameter_offset != columns {
899 return Err(TofMultiBankLatticeError::InternalInvariant);
900 }
901 let observed = bank
902 .input
903 .pattern
904 .observed_y
905 .as_ref()
906 .ok_or(TofLeBailError::MissingObservations)?;
907 for sample in 0..samples {
908 if bank
909 .input
910 .pattern
911 .mask
912 .as_ref()
913 .is_some_and(|mask| !mask[sample])
914 {
915 continue;
916 }
917 let weight = if options.use_uncertainty {
918 bank.input
919 .pattern
920 .uncertainty
921 .as_ref()
922 .map_or(1.0, |sigma| sigma[sample].recip())
923 } else {
924 1.0
925 };
926 selected_residual[row_offset] = (observed[sample] - calculation.y[sample]) * weight;
927 for parameter in 0..columns {
928 selected_jacobian[(row_offset, parameter)] = physical[(sample, parameter)] * weight;
929 }
930 row_offset += 1;
931 }
932 }
933 if row_offset != rows {
934 return Err(TofMultiBankLatticeError::InternalInvariant);
935 }
936 Ok((selected_jacobian, selected_residual))
937}
938
939pub(crate) fn lattice_changes(
940 input: &TofMultiBankLatticeInput,
941 before: &[f64],
942 after: &[f64],
943 scales: &[f64],
944) -> Result<Vec<TofLatticeParameterChange>, TofMultiBankLatticeError> {
945 if before.len() != after.len() || before.len() != scales.len() {
946 return Err(TofMultiBankLatticeError::InternalInvariant);
947 }
948 let mut changes = Vec::new();
949 let mut offset = 0;
950 for phase in &input.lattice_phases {
951 for name in phase.parameterization.parameter_names() {
952 if before[offset].to_bits() != after[offset].to_bits() {
953 changes.push(TofLatticeParameterChange {
954 phase_id: phase.phase_id.clone(),
955 parameter_name: name.clone(),
956 before: before[offset],
957 after: after[offset],
958 scaled_change: (after[offset] - before[offset]) / scales[offset],
959 });
960 }
961 offset += 1;
962 }
963 }
964 Ok(changes)
965}
966
967fn checkpoint_from_state(
968 input: &TofMultiBankLatticeInput,
969 states: &[AcceptedBankState],
970 cells: &[UnitCell],
971 history: &[TofMultiBankLatticeIterationRecord],
972) -> TofMultiBankLatticeCheckpoint {
973 TofMultiBankLatticeCheckpoint {
974 completed_iterations: history.len(),
975 banks: input
976 .multibank
977 .banks
978 .iter()
979 .zip(states)
980 .map(|(bank, state)| TofMultiBankCheckpointBank {
981 bank_id: bank.bank_id.clone(),
982 phases: state.phases.clone(),
983 background: state.background.clone(),
984 })
985 .collect(),
986 lattice_phases: lattice_states(input, cells),
987 history: history.to_vec(),
988 }
989}
990
991fn lattice_states(
992 input: &TofMultiBankLatticeInput,
993 cells: &[UnitCell],
994) -> Vec<TofSharedLatticeState> {
995 input
996 .lattice_phases
997 .iter()
998 .zip(cells)
999 .map(|(phase, cell)| TofSharedLatticeState {
1000 phase_id: phase.phase_id.clone(),
1001 cell: *cell,
1002 })
1003 .collect()
1004}
1005
1006fn validate_checkpoint_bank(
1007 saved: &TofMultiBankCheckpointBank,
1008 original: &crate::TofLeBailBank,
1009 input: &TofMultiBankLatticeInput,
1010 cells: &[TofSharedLatticeState],
1011) -> Result<(), TofMultiBankLatticeError> {
1012 if saved.bank_id != original.bank_id || saved.phases.len() != original.input.phases.len() {
1013 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
1014 "shared TOF lattice checkpoint bank identity or phase count changed",
1015 ));
1016 }
1017 for (saved_phase, original_phase) in saved.phases.iter().zip(&original.input.phases) {
1018 saved_phase.validate()?;
1019 if saved_phase.phase_id() != original_phase.phase_id()
1020 || saved_phase.name() != original_phase.name()
1021 || saved_phase.reflection_ids() != original_phase.reflection_ids()
1022 || saved_phase.hkl() != original_phase.hkl()
1023 || saved_phase.scale().to_bits() != original_phase.scale().to_bits()
1024 {
1025 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
1026 "shared TOF lattice checkpoint bank topology changed",
1027 ));
1028 }
1029 if let Some((model, cell)) = input
1030 .lattice_phases
1031 .iter()
1032 .zip(cells)
1033 .find(|(model, _)| model.phase_id() == saved_phase.phase_id())
1034 {
1035 let expected = tof_lattice_geometry(
1036 model.parameterization(),
1037 cell.cell,
1038 saved_phase.hkl(),
1039 original.input.instrument,
1040 )?;
1041 if saved_phase.d_spacing_angstrom() != expected.d_spacing_angstrom {
1042 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
1043 "shared TOF lattice checkpoint d-spacings disagree with its cell",
1044 ));
1045 }
1046 } else if saved_phase.d_spacing_angstrom() != original_phase.d_spacing_angstrom() {
1047 return Err(TofMultiBankLatticeError::InvalidCheckpoint(
1048 "fixed TOF phase d-spacings changed in a lattice checkpoint",
1049 ));
1050 }
1051 }
1052 match (&saved.background, &original.input.background) {
1053 (None, None) => Ok(()),
1054 (Some(saved), Some(original))
1055 if saved.background_id() == original.background_id()
1056 && saved
1057 .domain_us()
1058 .iter()
1059 .zip(original.domain_us())
1060 .all(|(saved, original)| saved.to_bits() == original.to_bits())
1061 && saved.coefficients().len() == original.coefficients().len() =>
1062 {
1063 Ok(())
1064 }
1065 _ => Err(TofMultiBankLatticeError::InvalidCheckpoint(
1066 "shared TOF lattice checkpoint background contract changed",
1067 )),
1068 }
1069}
1070
1071#[derive(Debug)]
1073pub enum TofMultiBankLatticeError {
1074 MultiBank(TofMultiBankError),
1076 Lattice(LatticeError),
1078 InvalidModel(&'static str),
1080 InvalidOptions,
1082 LinearSolve,
1084 InvalidCheckpoint(&'static str),
1086 AllocationOverflow,
1088 Runtime(RuntimeError),
1090 InternalInvariant,
1092}
1093
1094impl Display for TofMultiBankLatticeError {
1095 fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
1096 match self {
1097 Self::MultiBank(error) => Display::fmt(error, formatter),
1098 Self::Lattice(error) => Display::fmt(error, formatter),
1099 Self::InvalidModel(message) | Self::InvalidCheckpoint(message) => {
1100 formatter.write_str(message)
1101 }
1102 Self::InvalidOptions => formatter.write_str("invalid shared TOF lattice options"),
1103 Self::LinearSolve => formatter.write_str("shared TOF lattice linear solve failed"),
1104 Self::AllocationOverflow => {
1105 formatter.write_str("shared TOF lattice allocation overflow")
1106 }
1107 Self::Runtime(error) => Display::fmt(error, formatter),
1108 Self::InternalInvariant => {
1109 formatter.write_str("shared TOF lattice internal array invariant failed")
1110 }
1111 }
1112 }
1113}
1114
1115impl Error for TofMultiBankLatticeError {
1116 fn source(&self) -> Option<&(dyn Error + 'static)> {
1117 match self {
1118 Self::MultiBank(error) => Some(error),
1119 Self::Lattice(error) => Some(error),
1120 Self::Runtime(error) => Some(error),
1121 _ => None,
1122 }
1123 }
1124}
1125
1126impl From<TofMultiBankError> for TofMultiBankLatticeError {
1127 fn from(value: TofMultiBankError) -> Self {
1128 Self::MultiBank(value)
1129 }
1130}
1131
1132impl From<TofLeBailError> for TofMultiBankLatticeError {
1133 fn from(value: TofLeBailError) -> Self {
1134 Self::MultiBank(TofMultiBankError::LeBail(value))
1135 }
1136}
1137
1138impl From<LatticeError> for TofMultiBankLatticeError {
1139 fn from(value: LatticeError) -> Self {
1140 Self::Lattice(value)
1141 }
1142}
1143
1144impl From<RuntimeError> for TofMultiBankLatticeError {
1145 fn from(value: RuntimeError) -> Self {
1146 Self::Runtime(value)
1147 }
1148}