1use std::fmt;
11use std::sync::Arc;
12
13use nereids_core::constants::{EV_TO_JOULES, NEUTRON_MASS_KG, PIVOT_FLOOR, tof_to_energy};
14use nereids_endf::resonance::ResonanceData;
15use nereids_fitting::exact_count_model::ExactTwoArmRatioModel;
16use nereids_fitting::joint_poisson::{self, JointPoissonFitConfig, JointPoissonObjective};
17use nereids_fitting::lm::{self, FitModel, LmConfig, LmResult};
18use nereids_fitting::parameters::{FitParameter, ParameterSet};
19use nereids_fitting::poisson::{self, PoissonConfig};
20use nereids_fitting::transmission_model::{
21 EnergyScaleTransmissionModel, MultiplicativeBaselineModel, NormalizedTransmissionModel,
22 PrecomputedTransmissionModel, TransmissionFitModel,
23};
24use nereids_physics::counts_response::DetectorBinResponseMatrix;
25use nereids_physics::resolution::ResolutionFunction;
26use nereids_physics::transmission::{InstrumentParams, WorkingGridLayout, WorkingGridXs};
27
28use crate::error::PipelineError;
29
30#[derive(Debug, Clone)]
31pub struct PrecomputedXs {
32 pub sigma: Arc<Vec<Vec<f64>>>,
33 pub layout: Arc<WorkingGridLayout>,
34}
35
36impl From<WorkingGridXs> for PrecomputedXs {
37 fn from(working: WorkingGridXs) -> Self {
38 Self {
39 sigma: Arc::new(working.sigma),
40 layout: Arc::new(working.layout),
41 }
42 }
43}
44
45#[derive(Debug, Clone)]
59pub struct BackgroundConfig {
60 pub anorm_init: f64,
62 pub back_a_init: f64,
64 pub back_b_init: f64,
66 pub back_c_init: f64,
68 pub back_d_init: f64,
73 pub back_f_init: f64,
78 pub fit_anorm: bool,
80 pub fit_back_a: bool,
82 pub fit_back_b: bool,
84 pub fit_back_c: bool,
86 pub fit_back_d: bool,
88 pub fit_back_f: bool,
90}
91
92impl Default for BackgroundConfig {
93 fn default() -> Self {
94 Self {
95 anorm_init: 1.0,
96 back_a_init: 0.0,
97 back_b_init: 0.0,
98 back_c_init: 0.0,
99 back_d_init: 0.01,
100 back_f_init: 1.0,
101 fit_anorm: true,
102 fit_back_a: true,
103 fit_back_b: true,
104 fit_back_c: true,
105 fit_back_d: false,
106 fit_back_f: false,
107 }
108 }
109}
110
111#[derive(Debug, Clone, Copy)]
116struct BackgroundIndices {
117 anorm: usize,
118 back_a: usize,
119 back_b: usize,
120 back_c: usize,
121 back_d: Option<usize>,
122 back_f: Option<usize>,
123}
124
125#[derive(Debug, Clone)]
157pub struct MultiplicativeBaselineConfig {
158 pub b0_init: f64,
160 pub b1_init: f64,
162 pub b2_init: f64,
164 pub fit_b0: bool,
167 pub fit_b1: bool,
169 pub fit_b2: bool,
171 pub spatial_global: bool,
177 pub b0_bounds: (f64, f64),
179 pub b1_bounds: (f64, f64),
181 pub b2_bounds: (f64, f64),
183}
184
185impl Default for MultiplicativeBaselineConfig {
186 fn default() -> Self {
187 Self {
188 b0_init: 1.0,
189 b1_init: 0.0,
190 b2_init: 0.0,
191 fit_b0: true,
192 fit_b1: true,
193 fit_b2: true,
194 spatial_global: true,
195 b0_bounds: (0.9, 1.1),
196 b1_bounds: (-0.05, 0.05),
197 b2_bounds: (-0.05, 0.05),
198 }
199 }
200}
201
202#[derive(Debug, Clone, Copy)]
205struct BaselineIndices {
206 b0: usize,
207 b1: usize,
208 b2: usize,
209}
210
211#[derive(Debug, Clone)]
223pub enum InputData {
224 Transmission {
228 transmission: Vec<f64>,
230 uncertainty: Vec<f64>,
232 },
233 Counts {
239 sample_counts: Vec<f64>,
241 open_beam_counts: Vec<f64>,
243 },
244 CountsWithNuisance {
248 sample_counts: Vec<f64>,
250 flux: Vec<f64>,
252 background: Vec<f64>,
254 },
255}
256
257impl InputData {
258 pub fn n_energies(&self) -> usize {
260 match self {
261 Self::Transmission { transmission, .. } => transmission.len(),
262 Self::Counts { sample_counts, .. } => sample_counts.len(),
263 Self::CountsWithNuisance { sample_counts, .. } => sample_counts.len(),
264 }
265 }
266
267 pub fn is_counts(&self) -> bool {
269 matches!(self, Self::Counts { .. } | Self::CountsWithNuisance { .. })
270 }
271}
272
273#[derive(Debug, Clone, Default)]
278pub enum SolverConfig {
279 LevenbergMarquardt(LmConfig),
281 PoissonKL(PoissonConfig),
296 #[default]
299 Auto,
300}
301
302#[derive(Debug, Clone)]
331pub struct CountsBackgroundConfig {
332 pub alpha_1_init: f64,
338 pub alpha_2_init: f64,
341 pub fit_alpha_1: bool,
345 pub fit_alpha_2: bool,
347 pub c: f64,
356}
357
358impl Default for CountsBackgroundConfig {
359 fn default() -> Self {
360 Self {
361 alpha_1_init: 1.0,
362 alpha_2_init: 1.0,
363 fit_alpha_1: false,
364 fit_alpha_2: false,
365 c: 1.0,
366 }
367 }
368}
369
370#[derive(Debug, Clone)]
376pub struct ExactCountResponseConfig {
377 pub incident_fluence_weights: Vec<f64>,
383 pub detector_time_edges_us: Vec<f64>,
385 pub timing_offset_us: f64,
387}
388
389#[derive(Debug, Clone)]
396pub struct UnifiedFitConfig {
397 energies: Vec<f64>,
399 resonance_data: Vec<ResonanceData>,
400 isotope_names: Vec<String>,
401 temperature_k: f64,
402 resolution: Option<ResolutionFunction>,
403 initial_densities: Vec<f64>,
404 fit_temperature: bool,
405 compute_covariance: bool,
406 scale_by_chi2: bool,
411
412 solver: SolverConfig,
414
415 transmission_background: Option<BackgroundConfig>,
418 multiplicative_baseline: Option<MultiplicativeBaselineConfig>,
420 counts_background: Option<CountsBackgroundConfig>,
422 exact_count_response: Option<ExactCountResponseConfig>,
424
425 counts_enable_polish: Option<bool>,
432
433 precomputed_cross_sections: Option<PrecomputedXs>,
435 precomputed_base_xs: Option<Arc<Vec<Vec<f64>>>>,
436 precomputed_resolution_plan: Option<Arc<nereids_physics::resolution::ResolutionPlan>>,
445 precomputed_sparse_cubature_plan:
456 Option<Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>>,
457 precomputed_sparse_scalar_plan: Option<Arc<nereids_physics::surrogate::ScalarSurrogatePlan>>,
461
462 fit_energy_scale: bool,
467 t0_init_us: f64,
469 l_scale_init: f64,
471 flight_path_m: f64,
473 tzero_jacobian_method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,
479 energy_scale_seed_enabled: bool,
487
488 fit_energy_range: Option<(f64, f64)>,
508
509 pub(crate) density_indices: Option<Vec<usize>>,
513 pub(crate) density_ratios: Option<Vec<f64>>,
516 n_density_params: Option<usize>,
519
520 density_free: Option<Vec<bool>>,
537}
538
539impl UnifiedFitConfig {
540 pub fn new(
542 energies: Vec<f64>,
543 resonance_data: Vec<ResonanceData>,
544 isotope_names: Vec<String>,
545 temperature_k: f64,
546 resolution: Option<ResolutionFunction>,
547 initial_densities: Vec<f64>,
548 ) -> Result<Self, FitConfigError> {
549 if energies.is_empty() {
550 return Err(FitConfigError::EmptyEnergies);
551 }
552 if resonance_data.is_empty() {
553 return Err(FitConfigError::EmptyResonanceData);
554 }
555 if initial_densities.len() != resonance_data.len() {
556 return Err(FitConfigError::DensityCountMismatch {
557 densities: initial_densities.len(),
558 isotopes: resonance_data.len(),
559 });
560 }
561 if isotope_names.len() != resonance_data.len() {
562 return Err(FitConfigError::NameCountMismatch {
563 names: isotope_names.len(),
564 isotopes: resonance_data.len(),
565 });
566 }
567 if !temperature_k.is_finite() {
568 return Err(FitConfigError::NonFiniteTemperature(temperature_k));
569 }
570 if temperature_k < 0.0 {
571 return Err(FitConfigError::NegativeTemperature(temperature_k));
572 }
573 Ok(Self {
574 energies,
575 resonance_data,
576 isotope_names,
577 temperature_k,
578 resolution,
579 initial_densities,
580 fit_temperature: false,
581 compute_covariance: true,
582 scale_by_chi2: false,
583 solver: SolverConfig::Auto,
584 transmission_background: None,
585 multiplicative_baseline: None,
586 counts_background: None,
587 exact_count_response: None,
588 counts_enable_polish: None,
589 precomputed_cross_sections: None,
590 precomputed_base_xs: None,
591 precomputed_resolution_plan: None,
592 precomputed_sparse_cubature_plan: None,
593 precomputed_sparse_scalar_plan: None,
594 fit_energy_scale: false,
595 t0_init_us: 0.0,
596 l_scale_init: 1.0,
597 flight_path_m: 25.0,
598 density_indices: None,
599 density_ratios: None,
600 n_density_params: None,
601 density_free: None,
602 tzero_jacobian_method: None,
603 energy_scale_seed_enabled: true,
604 fit_energy_range: None,
605 })
606 }
607
608 #[must_use]
614 pub fn with_tzero_jacobian_method(
615 mut self,
616 method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,
617 ) -> Self {
618 self.tzero_jacobian_method = method;
619 self
620 }
621
622 #[must_use]
629 pub fn with_energy_scale_seed(mut self, enabled: bool) -> Self {
630 self.energy_scale_seed_enabled = enabled;
631 self
632 }
633
634 #[must_use]
637 pub fn with_solver(mut self, solver: SolverConfig) -> Self {
638 self.solver = solver;
639 self
640 }
641
642 #[must_use]
643 pub fn with_fit_temperature(mut self, v: bool) -> Self {
644 self.fit_temperature = v;
645 self
646 }
647
648 #[must_use]
649 pub fn with_compute_covariance(mut self, v: bool) -> Self {
650 self.compute_covariance = v;
651 self
652 }
653
654 #[must_use]
671 pub fn with_scale_by_chi2(mut self, v: bool) -> Self {
672 self.scale_by_chi2 = v;
673 self
674 }
675
676 #[must_use]
682 pub fn with_energy_scale(
683 mut self,
684 t0_init_us: f64,
685 l_scale_init: f64,
686 flight_path_m: f64,
687 ) -> Self {
688 self.fit_energy_scale = true;
689 self.t0_init_us = t0_init_us;
690 self.l_scale_init = l_scale_init;
691 self.flight_path_m = flight_path_m;
692 self
693 }
694
695 pub fn with_fit_energy_range(
708 mut self,
709 range: Option<(f64, f64)>,
710 ) -> Result<Self, FitConfigError> {
711 if let Some((lo, hi)) = range {
712 if !lo.is_finite() || !hi.is_finite() {
713 return Err(FitConfigError::InvalidFitEnergyRange(
714 "fit_energy_range bounds must be finite",
715 ));
716 }
717 if lo >= hi {
718 return Err(FitConfigError::InvalidFitEnergyRange(
719 "fit_energy_range min must be strictly less than max",
720 ));
721 }
722 }
723 self.fit_energy_range = range;
724 Ok(self)
725 }
726
727 #[must_use]
728 pub fn with_transmission_background(mut self, bg: BackgroundConfig) -> Self {
729 self.transmission_background = Some(bg);
730 self
731 }
732
733 #[must_use]
740 pub fn with_multiplicative_baseline(mut self, bl: MultiplicativeBaselineConfig) -> Self {
741 self.multiplicative_baseline = Some(bl);
742 self
743 }
744
745 #[must_use]
746 pub fn with_counts_background(mut self, bg: CountsBackgroundConfig) -> Self {
747 self.counts_background = Some(bg);
748 self
749 }
750
751 #[must_use]
753 pub fn with_exact_count_response(mut self, response: ExactCountResponseConfig) -> Self {
754 self.exact_count_response = Some(response);
755 self
756 }
757
758 #[must_use]
763 pub fn with_counts_enable_polish(mut self, v: Option<bool>) -> Self {
764 self.counts_enable_polish = v;
765 self
766 }
767
768 #[must_use]
769 pub fn with_precomputed_cross_sections(mut self, xs: PrecomputedXs) -> Self {
770 self.precomputed_cross_sections = Some(xs);
771 self.precomputed_sparse_cubature_plan = None;
772 self.precomputed_sparse_scalar_plan = None;
773 self
774 }
775
776 #[must_use]
777 pub fn with_precomputed_base_xs(mut self, xs: Arc<Vec<Vec<f64>>>) -> Self {
778 self.precomputed_base_xs = Some(xs);
779 self.precomputed_sparse_cubature_plan = None;
783 self.precomputed_sparse_scalar_plan = None;
784 self
785 }
786
787 #[must_use]
795 pub fn with_precomputed_resolution_plan(
796 mut self,
797 plan: Arc<nereids_physics::resolution::ResolutionPlan>,
798 ) -> Self {
799 self.precomputed_resolution_plan = Some(plan);
800 self
801 }
802
803 #[must_use]
821 pub fn with_precomputed_sparse_cubature_plan(
822 mut self,
823 plan: Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>,
824 ) -> Self {
825 self.precomputed_sparse_cubature_plan = Some(plan);
826 self
827 }
828
829 #[must_use]
835 pub fn with_precomputed_sparse_scalar_plan(
836 mut self,
837 plan: Arc<nereids_physics::surrogate::ScalarSurrogatePlan>,
838 ) -> Self {
839 self.precomputed_sparse_scalar_plan = Some(plan);
840 self
841 }
842
843 pub fn with_groups(
860 mut self,
861 groups: &[(&nereids_core::types::IsotopeGroup, &[ResonanceData])],
862 initial_densities: Vec<f64>,
863 ) -> Result<Self, FitConfigError> {
864 if self.density_free.is_some() {
865 return Err(FitConfigError::DensityFreezeBeforeGroups);
866 }
867 if groups.is_empty() {
868 return Err(FitConfigError::EmptyResonanceData);
869 }
870 if initial_densities.len() != groups.len() {
871 return Err(FitConfigError::DensityCountMismatch {
872 densities: initial_densities.len(),
873 isotopes: groups.len(),
874 });
875 }
876 let mut all_resonance_data = Vec::new();
877 let mut all_indices = Vec::new();
878 let mut all_ratios = Vec::new();
879 let mut names = Vec::new();
880 for (g_idx, (group, rd_list)) in groups.iter().enumerate() {
881 if rd_list.len() != group.n_members() {
882 return Err(FitConfigError::GroupMemberCountMismatch {
883 group_name: group.name().to_string(),
884 rd_count: rd_list.len(),
885 member_count: group.n_members(),
886 });
887 }
888 names.push(group.name().to_string());
889 for ((isotope, ratio), rd) in group.members().iter().zip(rd_list.iter()) {
890 if rd.isotope != *isotope {
892 return Err(FitConfigError::GroupMemberIsotopeMismatch {
893 group_name: group.name().to_string(),
894 expected_z: isotope.z(),
895 expected_a: isotope.a(),
896 got_z: rd.isotope.z(),
897 got_a: rd.isotope.a(),
898 });
899 }
900 all_resonance_data.push(rd.clone());
901 all_indices.push(g_idx);
902 all_ratios.push(*ratio);
903 }
904 }
905 self.resonance_data = all_resonance_data;
906 self.isotope_names = names;
907 self.initial_densities = initial_densities;
908 self.n_density_params = Some(groups.len());
909 self.density_indices = Some(all_indices);
910 self.density_ratios = Some(all_ratios);
911 self.precomputed_cross_sections = None;
917 self.precomputed_base_xs = None;
918 self.precomputed_sparse_cubature_plan = None;
924 self.precomputed_sparse_scalar_plan = None;
925 Ok(self)
926 }
927
928 pub fn precomputed_sparse_cubature_plan(
935 &self,
936 ) -> Option<&Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>> {
937 self.precomputed_sparse_cubature_plan.as_ref()
938 }
939
940 pub fn precomputed_sparse_scalar_plan(
942 &self,
943 ) -> Option<&Arc<nereids_physics::surrogate::ScalarSurrogatePlan>> {
944 self.precomputed_sparse_scalar_plan.as_ref()
945 }
946
947 pub fn energies(&self) -> &[f64] {
948 &self.energies
949 }
950 pub fn resonance_data(&self) -> &[ResonanceData] {
951 &self.resonance_data
952 }
953 pub fn isotope_names(&self) -> &[String] {
954 &self.isotope_names
955 }
956 pub fn temperature_k(&self) -> f64 {
957 self.temperature_k
958 }
959 pub fn resolution(&self) -> Option<&ResolutionFunction> {
960 self.resolution.as_ref()
961 }
962 pub fn initial_densities(&self) -> &[f64] {
963 &self.initial_densities
964 }
965 pub fn solver(&self) -> &SolverConfig {
966 &self.solver
967 }
968 pub fn fit_temperature(&self) -> bool {
969 self.fit_temperature
970 }
971 pub fn transmission_background(&self) -> Option<&BackgroundConfig> {
972 self.transmission_background.as_ref()
973 }
974 pub fn multiplicative_baseline(&self) -> Option<&MultiplicativeBaselineConfig> {
976 self.multiplicative_baseline.as_ref()
977 }
978 pub fn counts_background(&self) -> Option<&CountsBackgroundConfig> {
979 self.counts_background.as_ref()
980 }
981 pub fn exact_count_response(&self) -> Option<&ExactCountResponseConfig> {
982 self.exact_count_response.as_ref()
983 }
984 pub fn counts_enable_polish(&self) -> Option<bool> {
986 self.counts_enable_polish
987 }
988 pub fn fit_energy_scale(&self) -> bool {
991 self.fit_energy_scale
992 }
993 pub fn scale_by_chi2(&self) -> bool {
996 self.scale_by_chi2
997 }
998 pub fn flight_path_m(&self) -> f64 {
1000 self.flight_path_m
1001 }
1002 pub fn fit_energy_range(&self) -> Option<(f64, f64)> {
1006 self.fit_energy_range
1007 }
1008 pub fn baseline_reference_energy(&self) -> f64 {
1016 let mask = nereids_fitting::active_mask::build_active_mask(
1017 self.energies(),
1018 self.fit_energy_range(),
1019 );
1020 nereids_fitting::transmission_model::baseline_reference_energy_active(
1021 self.energies(),
1022 mask.as_deref(),
1023 )
1024 }
1025 pub fn precomputed_cross_sections(&self) -> Option<&PrecomputedXs> {
1026 self.precomputed_cross_sections.as_ref()
1027 }
1028 pub fn n_density_params(&self) -> usize {
1030 self.n_density_params.unwrap_or(self.resonance_data.len())
1031 }
1032
1033 #[must_use]
1046 pub fn with_fix_densities(mut self, fix: bool) -> Self {
1047 self.density_free = if fix {
1048 Some(vec![false; self.n_density_params()])
1049 } else {
1050 None
1051 };
1052 self
1053 }
1054
1055 pub fn with_density_free(mut self, free: Vec<bool>) -> Result<Self, FitConfigError> {
1069 let n = self.n_density_params();
1070 if free.len() != n {
1071 return Err(FitConfigError::DensityCountMismatch {
1072 densities: free.len(),
1073 isotopes: n,
1074 });
1075 }
1076 self.density_free = if free.iter().all(|&f| f) {
1079 None
1080 } else {
1081 Some(free)
1082 };
1083 Ok(self)
1084 }
1085
1086 fn density_is_fixed(&self, i: usize) -> bool {
1089 self.density_free
1090 .as_ref()
1091 .is_some_and(|mask| !mask.get(i).copied().unwrap_or(true))
1092 }
1093
1094 fn n_free_density_params(&self) -> usize {
1097 (0..self.n_density_params())
1098 .filter(|&i| !self.density_is_fixed(i))
1099 .count()
1100 }
1101
1102 pub(crate) fn effective_solver(&self, input: &InputData) -> SolverConfig {
1104 match &self.solver {
1105 SolverConfig::Auto => {
1106 if input.is_counts() {
1107 SolverConfig::PoissonKL(PoissonConfig::default())
1108 } else {
1109 SolverConfig::LevenbergMarquardt(LmConfig::default())
1110 }
1111 }
1112 other => other.clone(),
1113 }
1114 }
1115}
1116
1117pub(crate) fn validate_counts_resolution_route(
1128 is_counts: bool,
1129 observed_bin_count: usize,
1130 config: &UnifiedFitConfig,
1131) -> Result<(), PipelineError> {
1132 if !is_counts && config.exact_count_response().is_some() {
1133 return Err(PipelineError::InvalidParameter(
1134 "exact_count_response is a raw-count model and cannot be attached to \
1135 normalized transmission"
1136 .into(),
1137 ));
1138 }
1139 if is_counts && let Some(resolution) = config.resolution() {
1140 if matches!(resolution, ResolutionFunction::Gaussian(_)) {
1144 return Err(PipelineError::InvalidParameter(
1145 "counts input with instrument resolution requires the exact \
1146 separate-arm model R[Phi] and R[Phi*T], and Gaussian energy \
1147 broadening cannot provide it: resolved raw-count fitting needs \
1148 a detector-time response (TabulatedResolution or \
1149 IkedaCarpenter); the R[T] shortcut is not used"
1150 .into(),
1151 ));
1152 }
1153 let exact = config.exact_count_response().ok_or_else(|| {
1154 PipelineError::InvalidParameter(
1155 "counts input with instrument resolution requires the exact \
1156 separate-arm model R[Phi] and R[Phi*T]: provide incident fluence \
1157 weights and detector-time bin edges through exact_count_response; \
1158 the R[T] shortcut is not used"
1159 .into(),
1160 )
1161 })?;
1162 if exact.incident_fluence_weights.len() != config.energies().len() {
1163 return Err(PipelineError::ShapeMismatch(format!(
1164 "exact_count_response incident_fluence_weights length {} must match \
1165 the true-energy grid length {}",
1166 exact.incident_fluence_weights.len(),
1167 config.energies().len()
1168 )));
1169 }
1170 for (index, &fluence) in exact.incident_fluence_weights.iter().enumerate() {
1177 if !fluence.is_finite() || fluence < 0.0 {
1178 return Err(PipelineError::InvalidParameter(format!(
1179 "exact_count_response incident_fluence_weights[{index}] must be \
1180 finite and >= 0, got {fluence}"
1181 )));
1182 }
1183 }
1184 if !exact.incident_fluence_weights.iter().any(|&f| f > 0.0) {
1185 return Err(PipelineError::InvalidParameter(
1186 "exact_count_response incident_fluence_weights must contain at least \
1187 one positive value; an all-zero incident source cannot produce the \
1188 observed counts"
1189 .into(),
1190 ));
1191 }
1192 if exact.detector_time_edges_us.len() != observed_bin_count + 1 {
1193 return Err(PipelineError::ShapeMismatch(format!(
1194 "exact_count_response detector_time_edges_us length {} must be one \
1195 greater than the observed count-bin length {}",
1196 exact.detector_time_edges_us.len(),
1197 observed_bin_count
1198 )));
1199 }
1200 if !exact.timing_offset_us.is_finite() {
1201 return Err(PipelineError::InvalidParameter(format!(
1202 "exact_count_response timing_offset_us must be finite, got {}",
1203 exact.timing_offset_us
1204 )));
1205 }
1206 if config.fit_energy_scale() {
1207 return Err(PipelineError::InvalidParameter(
1208 "energy-scale fitting is not yet connected to the exact two-arm \
1209 detector response; fitting cross-section energy while holding the \
1210 response clock fixed would be physically inconsistent"
1211 .into(),
1212 ));
1213 }
1214 if config.fit_energy_range().is_some() {
1215 return Err(PipelineError::InvalidParameter(
1216 "fit_energy_range is not yet supported by exact detector-time count \
1217 response because true-energy points and detector-time bins are \
1218 different axes"
1219 .into(),
1220 ));
1221 }
1222 } else if is_counts && config.exact_count_response().is_some() {
1223 return Err(PipelineError::InvalidParameter(
1224 "exact_count_response requires an instrument resolution model".into(),
1225 ));
1226 }
1227 Ok(())
1228}
1229
1230pub fn fit_spectrum_typed(
1243 input: &InputData,
1244 config: &UnifiedFitConfig,
1245) -> Result<SpectrumFitResult, PipelineError> {
1246 validate_precomputed_cross_sections(config)?;
1247 fit_spectrum_validated(input, config)
1248}
1249
1250pub(crate) fn fit_spectrum_validated(
1251 input: &InputData,
1252 config: &UnifiedFitConfig,
1253) -> Result<SpectrumFitResult, PipelineError> {
1254 let n_e = config.energies().len();
1255
1256 if config.fit_temperature && config.temperature_k < 1.0 {
1258 return Err(PipelineError::InvalidParameter(format!(
1259 "temperature must be >= 1.0 K when fit_temperature is true, got {}",
1260 config.temperature_k,
1261 )));
1262 }
1263
1264 if count_free_params(config) == 0 {
1270 return Err(PipelineError::InvalidParameter(
1271 "no free parameters to fit: all densities are frozen and no other \
1272 parameter is free — free at least one density (with_density_free) \
1273 or enable fit_temperature / energy-scale / background"
1274 .to_string(),
1275 ));
1276 }
1277
1278 let exact_resolved_counts = input.is_counts() && config.exact_count_response().is_some();
1282 if input.n_energies() != n_e && !exact_resolved_counts {
1283 return Err(PipelineError::ShapeMismatch(format!(
1284 "input data has {} energy bins but config.energies has {}",
1285 input.n_energies(),
1286 n_e,
1287 )));
1288 }
1289
1290 match input {
1292 InputData::Transmission {
1293 transmission,
1294 uncertainty,
1295 } => {
1296 if uncertainty.len() != transmission.len() {
1297 return Err(PipelineError::ShapeMismatch(format!(
1298 "uncertainty length {} != transmission length {}",
1299 uncertainty.len(),
1300 transmission.len(),
1301 )));
1302 }
1303 }
1304 InputData::Counts {
1305 sample_counts,
1306 open_beam_counts,
1307 } => {
1308 if open_beam_counts.len() != sample_counts.len() {
1309 return Err(PipelineError::ShapeMismatch(format!(
1310 "open_beam_counts length {} != sample_counts length {}",
1311 open_beam_counts.len(),
1312 sample_counts.len(),
1313 )));
1314 }
1315 }
1316 InputData::CountsWithNuisance {
1317 sample_counts,
1318 flux,
1319 background,
1320 } => {
1321 if flux.len() != sample_counts.len() {
1322 return Err(PipelineError::ShapeMismatch(format!(
1323 "flux length {} != sample_counts length {}",
1324 flux.len(),
1325 sample_counts.len(),
1326 )));
1327 }
1328 if background.len() != sample_counts.len() {
1329 return Err(PipelineError::ShapeMismatch(format!(
1330 "background length {} != sample_counts length {}",
1331 background.len(),
1332 sample_counts.len(),
1333 )));
1334 }
1335 }
1336 }
1337
1338 if matches!(input, InputData::Transmission { .. }) && config.counts_background().is_some() {
1339 return Err(PipelineError::InvalidParameter(
1340 "counts background configuration cannot be used with transmission data; \
1341 use SAMMY transmission_background or multiplicative_baseline for a \
1342 transmission fit, or supply separate open/sample counts"
1343 .into(),
1344 ));
1345 }
1346
1347 validate_counts_resolution_route(input.is_counts(), input.n_energies(), config)?;
1348
1349 let effective_solver = config.effective_solver(input);
1350
1351 match (input, &effective_solver) {
1352 (
1354 InputData::Transmission {
1355 transmission,
1356 uncertainty,
1357 },
1358 SolverConfig::LevenbergMarquardt(lm_cfg),
1359 ) => fit_transmission_lm(transmission, uncertainty, config, lm_cfg),
1360
1361 (InputData::Transmission { .. }, SolverConfig::PoissonKL(_)) => {
1363 Err(PipelineError::InvalidParameter(
1364 "normalized transmission cannot use the Poisson/KL count objective: \
1365 fractional transmission is not Poisson count data and the supplied \
1366 uncertainty would be ignored; use the LM least-squares transmission \
1367 engine, or supply separate open/sample counts"
1368 .into(),
1369 ))
1370 }
1371
1372 (
1382 InputData::Counts {
1383 sample_counts,
1384 open_beam_counts,
1385 },
1386 SolverConfig::PoissonKL(poisson_cfg),
1387 ) => {
1388 let bg = vec![0.0f64; sample_counts.len()];
1389 fit_counts_joint_poisson(
1390 sample_counts,
1391 open_beam_counts,
1392 &bg,
1393 config,
1394 &poisson_to_joint_poisson_config(poisson_cfg, config),
1395 )
1396 }
1397
1398 (
1400 InputData::CountsWithNuisance {
1401 sample_counts,
1402 flux,
1403 background,
1404 },
1405 SolverConfig::PoissonKL(poisson_cfg),
1406 ) => fit_counts_joint_poisson(
1407 sample_counts,
1408 flux,
1409 background,
1410 config,
1411 &poisson_to_joint_poisson_config(poisson_cfg, config),
1412 ),
1413
1414 (InputData::Counts { .. }, SolverConfig::LevenbergMarquardt(_)) => {
1416 Err(PipelineError::InvalidParameter(
1417 "separate open/sample counts cannot use the LM least-squares \
1418 transmission engine: silently dividing the arms loses open-beam \
1419 uncertainty and count statistics; use the Poisson/KL count engine"
1420 .into(),
1421 ))
1422 }
1423
1424 (InputData::CountsWithNuisance { .. }, SolverConfig::LevenbergMarquardt(_)) => {
1426 Err(PipelineError::InvalidParameter(
1427 "CountsWithNuisance requires a counts-domain solver (LM cannot use nuisance parameters)"
1428 .into(),
1429 ))
1430 }
1431
1432 (_, SolverConfig::Auto) => unreachable!("Auto should be resolved before dispatch"),
1434 }
1435}
1436
1437fn poisson_to_joint_poisson_config(
1443 poisson_cfg: &PoissonConfig,
1444 config: &UnifiedFitConfig,
1445) -> JointPoissonFitConfig {
1446 let mut jp_cfg = JointPoissonFitConfig {
1447 max_iter: poisson_cfg.max_iter,
1448 scale_by_chi2: config.scale_by_chi2,
1449 ..Default::default()
1450 };
1451 if let Some(override_val) = config.counts_enable_polish() {
1452 jp_cfg.enable_polish = override_val;
1453 }
1454 jp_cfg
1455}
1456
1457fn fit_transmission_lm(
1459 measured_t: &[f64],
1460 sigma: &[f64],
1461 config: &UnifiedFitConfig,
1462 lm_config: &LmConfig,
1463) -> Result<SpectrumFitResult, PipelineError> {
1464 let n_density_params = config.n_density_params();
1465
1466 let mut param_vec = build_density_params(config);
1468
1469 let temperature_index = append_temperature_param(&mut param_vec, config);
1470 let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);
1471
1472 seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, measured_t, config);
1476
1477 if let Some(bg) = config.transmission_background.as_ref() {
1480 validate_transmission_background(bg)?;
1481 }
1482 let bg_indices = config
1483 .transmission_background
1484 .as_ref()
1485 .map(|bg| append_background_params(&mut param_vec, bg));
1486
1487 validate_multiplicative_baseline(config)?;
1490 let bl_indices = config
1491 .multiplicative_baseline
1492 .as_ref()
1493 .map(|bl| append_multiplicative_baseline_params(&mut param_vec, bl));
1494
1495 let mut params = ParameterSet::new(param_vec);
1496 let mut lm_cfg = lm_config.clone();
1497 lm_cfg.compute_covariance = config.compute_covariance;
1498
1499 let model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1501 build_energy_scale_transmission_model(config, t0_idx, ls_idx, temperature_index)?
1502 } else {
1503 build_transmission_model(config, n_density_params, temperature_index)?
1504 };
1505
1506 let active_mask = nereids_fitting::active_mask::build_active_mask(
1510 config.energies(),
1511 config.fit_energy_range(),
1512 );
1513
1514 if let Some((e_min, e_max)) = config.fit_energy_range() {
1525 let n_active = nereids_fitting::active_mask::active_count(
1526 active_mask.as_deref(),
1527 config.energies().len(),
1528 );
1529 let required = required_active_bins(config);
1530 if n_active < required {
1531 return Err(PipelineError::InvalidParameter(format!(
1532 "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
1533 on the configured energy grid; at least {required} active bin(s) are \
1534 required for LM transmission fitting with {n_free} free parameter(s) \
1535 (underdetermined when n_active < n_free)",
1536 n_free = count_free_params(config),
1537 )));
1538 }
1539 }
1540
1541 let mut stacked: Box<dyn FitModel> = model;
1547 if let Some(bi) = bg_indices {
1548 stacked = if let (Some(di), Some(fi)) = (bi.back_d, bi.back_f) {
1549 Box::new(NormalizedTransmissionModel::new_with_exponential(
1550 stacked,
1551 config.energies(),
1552 bi.anorm,
1553 bi.back_a,
1554 bi.back_b,
1555 bi.back_c,
1556 di,
1557 fi,
1558 ))
1559 } else {
1560 Box::new(NormalizedTransmissionModel::new(
1561 stacked,
1562 config.energies(),
1563 bi.anorm,
1564 bi.back_a,
1565 bi.back_b,
1566 bi.back_c,
1567 ))
1568 };
1569 }
1570 if let Some(bli) = bl_indices {
1571 stacked = Box::new(
1572 MultiplicativeBaselineModel::new(
1573 stacked,
1574 config.energies(),
1575 nereids_fitting::transmission_model::baseline_reference_energy_active(
1576 config.energies(),
1577 active_mask.as_deref(),
1578 ),
1579 bli.b0,
1580 bli.b1,
1581 bli.b2,
1582 )
1583 .with_active_mask(active_mask.as_deref()),
1585 );
1586 }
1587 let result = lm::levenberg_marquardt_with_mask(
1588 &*stacked,
1589 measured_t,
1590 sigma,
1591 &mut params,
1592 &lm_cfg,
1593 active_mask.as_deref(),
1594 )?;
1595
1596 let free_indices = params.free_indices();
1597 let mut sr = extract_result(
1598 config,
1599 &result,
1600 n_density_params,
1601 &free_indices,
1602 bg_indices,
1603 bl_indices,
1604 )?;
1605
1606 if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1608 sr.t0_us = Some(result.params[t0_idx]);
1609 sr.l_scale = Some(result.params[ls_idx]);
1610 sr.energy_scale_flight_path_m = Some(config.flight_path_m);
1611 }
1612
1613 Ok(sr)
1614}
1615
1616fn fit_counts_joint_poisson(
1642 sample_counts: &[f64],
1643 flux: &[f64],
1644 detector_background: &[f64],
1645 config: &UnifiedFitConfig,
1646 jp_cfg: &JointPoissonFitConfig,
1647) -> Result<SpectrumFitResult, PipelineError> {
1648 if let Some(bg) = config.counts_background()
1650 && (bg.fit_alpha_1 || bg.fit_alpha_2)
1651 {
1652 return Err(PipelineError::InvalidParameter(
1653 "joint-Poisson solver does not support fit_alpha_1/fit_alpha_2: \
1654 the profile lambda-hat absorbs the global flux scale (alpha_1 redundant); \
1655 alpha_2 / B_det wiring is not yet implemented."
1656 .into(),
1657 ));
1658 }
1659 if detector_background
1660 .iter()
1661 .any(|&v| !v.is_finite() || v < 0.0)
1662 {
1663 return Err(PipelineError::InvalidParameter(
1664 "detector_background must be finite and non-negative expected counts".into(),
1665 ));
1666 }
1667
1668 if let Some(bg) = config.transmission_background.as_ref() {
1670 if bg.fit_back_d || bg.fit_back_f {
1671 return Err(PipelineError::InvalidParameter(
1672 "joint-Poisson solver does not support the BackD/BackF exponential \
1673 tail (support is deferred)."
1674 .into(),
1675 ));
1676 }
1677 if (bg.fit_back_b || bg.fit_back_c) && !bg.fit_back_a {
1678 return Err(PipelineError::InvalidParameter(
1679 "joint-Poisson transmission_background: B_A (fit_back_a) must be \
1680 enabled whenever any of B_B / B_C is enabled (A_n alone cannot \
1681 absorb a constant offset — benchmarked at −23% density bias)."
1682 .into(),
1683 ));
1684 }
1685 }
1686
1687 let c = config.counts_background().map(|b| b.c).unwrap_or(1.0);
1688 if !(c.is_finite() && c > 0.0) {
1689 return Err(PipelineError::InvalidParameter(format!(
1690 "joint-Poisson solver requires finite c > 0 in CountsBackgroundConfig, got {c}",
1691 )));
1692 }
1693
1694 let n_density_params = config.n_density_params();
1699 let mut param_vec = build_density_params(config);
1700
1701 let temperature_index = append_temperature_param(&mut param_vec, config);
1702 let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);
1703
1704 if energy_scale_indices.is_some() {
1709 let t_proxy: Vec<f64> = sample_counts
1710 .iter()
1711 .zip(flux.iter())
1712 .zip(detector_background.iter())
1713 .map(|((&s, &f), &b)| {
1714 let den = f - b;
1715 if den > 0.0 {
1716 ((s - b).max(0.0) / den).min(2.0)
1717 } else {
1718 1.0
1719 }
1720 })
1721 .collect();
1722 seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, &t_proxy, config);
1723 }
1724
1725 let bg_indices = config
1729 .transmission_background
1730 .as_ref()
1731 .map(|bg| append_background_params(&mut param_vec, bg));
1732
1733 validate_multiplicative_baseline(config)?;
1738 let bl_indices = config
1739 .multiplicative_baseline
1740 .as_ref()
1741 .map(|bl| append_multiplicative_baseline_params(&mut param_vec, bl));
1742
1743 let mut params = ParameterSet::new(param_vec);
1744
1745 let mut true_model_config = config.clone();
1750 if config.exact_count_response().is_some() {
1751 true_model_config.resolution = None;
1752 true_model_config.precomputed_resolution_plan = None;
1753 true_model_config.precomputed_sparse_cubature_plan = None;
1754 true_model_config.precomputed_sparse_scalar_plan = None;
1755 }
1756 let t_model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1757 build_energy_scale_transmission_model(
1758 &true_model_config,
1759 t0_idx,
1760 ls_idx,
1761 temperature_index,
1762 )?
1763 } else {
1764 build_transmission_model(&true_model_config, n_density_params, temperature_index)?
1765 };
1766
1767 let mut active_mask = nereids_fitting::active_mask::build_active_mask(
1778 config.energies(),
1779 config.fit_energy_range(),
1780 );
1781
1782 if let Some((e_min, e_max)) = config.fit_energy_range() {
1792 let n_active = nereids_fitting::active_mask::active_count(
1793 active_mask.as_deref(),
1794 config.energies().len(),
1795 );
1796 let required = required_active_bins(config);
1797 if n_active < required {
1798 return Err(PipelineError::InvalidParameter(format!(
1799 "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
1800 on the configured energy grid; at least {required} active bin(s) are \
1801 required for joint-Poisson fitting with {n_free} free parameter(s) \
1802 (underdetermined when n_active < n_free)",
1803 n_free = count_free_params(config),
1804 )));
1805 }
1806 }
1807
1808 let mut stacked: Box<dyn FitModel> = t_model;
1814 if let Some(exact) = config.exact_count_response() {
1815 if let Some(bli) = bl_indices {
1816 stacked = Box::new(MultiplicativeBaselineModel::new(
1817 stacked,
1818 config.energies(),
1819 nereids_fitting::transmission_model::baseline_reference_energy_active(
1820 config.energies(),
1821 None,
1822 ),
1823 bli.b0,
1824 bli.b1,
1825 bli.b2,
1826 ));
1827 }
1828 let response = config.resolution().expect(
1829 "validate_counts_resolution_route requires resolution with exact_count_response",
1830 );
1831 let matrix = DetectorBinResponseMatrix::new(
1832 config.energies(),
1833 &exact.detector_time_edges_us,
1834 exact.timing_offset_us,
1835 response,
1836 )
1837 .map_err(|error| {
1838 PipelineError::InvalidParameter(format!(
1839 "failed to build exact detector-bin response: {error}"
1840 ))
1841 })?;
1842 let exact_model =
1848 ExactTwoArmRatioModel::new(stacked, matrix, &exact.incident_fluence_weights).map_err(
1849 |error| {
1850 PipelineError::InvalidParameter(format!(
1851 "failed to build the exact two-arm ratio model: {error}"
1852 ))
1853 },
1854 )?;
1855 let mut dead_occupied = Vec::new();
1867 for (bin, ((&observed_open, &observed_sample), &predicted_open)) in flux
1868 .iter()
1869 .zip(sample_counts)
1870 .zip(exact_model.open_expectation())
1871 .enumerate()
1872 {
1873 if observed_open + observed_sample > 0.0 && predicted_open <= PIVOT_FLOOR {
1874 if detector_background[bin] <= 0.0 {
1878 return Err(PipelineError::InvalidParameter(format!(
1879 "exact incident source has zero detector response in occupied \
1880 detector bin {bin}; the supplied source/response cannot explain \
1881 the observed counts"
1882 )));
1883 }
1884 dead_occupied.push(bin);
1885 }
1886 }
1887 if !dead_occupied.is_empty() {
1888 let mask = active_mask.get_or_insert_with(|| vec![true; flux.len()]);
1889 for &bin in &dead_occupied {
1890 mask[bin] = false;
1891 }
1892 let remaining = mask.iter().filter(|&&active| active).count();
1893 let required = required_active_bins(config);
1894 if remaining < required {
1895 return Err(PipelineError::InvalidParameter(format!(
1896 "{} detector bin(s) hold only background and cannot be fitted, \
1897 leaving {remaining} active bin(s); at least {required} are required",
1898 dead_occupied.len()
1899 )));
1900 }
1901 }
1902 stacked = Box::new(exact_model);
1903
1904 if let Some(bi) = bg_indices {
1905 let flight_path_m = response.flight_path_m();
1920 let raw_energies: Vec<Option<f64>> = exact
1921 .detector_time_edges_us
1922 .windows(2)
1923 .map(|edges| {
1924 let corrected_tof_us = 0.5 * (edges[0] + edges[1]) - exact.timing_offset_us;
1925 let energy = tof_to_energy(corrected_tof_us, flight_path_m);
1926 (energy.is_finite() && energy > 0.0).then_some(energy)
1927 })
1928 .collect();
1929 for (bin, energy) in raw_energies.iter().enumerate() {
1930 let occupied = flux[bin] + sample_counts[bin] > 0.0;
1931 if energy.is_none() && occupied {
1932 return Err(PipelineError::InvalidParameter(format!(
1933 "detector-time bin {bin} carries observed counts but has \
1934 non-positive time after the fixed timing offset, so the \
1935 SAMMY apparent-transmission background has no physical \
1936 energy there; check timing_offset_us against the \
1937 detector-time axis"
1938 )));
1939 }
1940 }
1941 let first_valid = raw_energies
1946 .iter()
1947 .flatten()
1948 .copied()
1949 .next()
1950 .ok_or_else(|| {
1951 PipelineError::InvalidParameter(
1952 "no detector-time bin has a positive time after the fixed timing \
1953 offset; the entire acquisition window precedes timing_offset_us"
1954 .to_string(),
1955 )
1956 })?;
1957 let mut last_valid = first_valid;
1958 let detector_energies: Vec<f64> = raw_energies
1959 .into_iter()
1960 .map(|energy| {
1961 if let Some(value) = energy {
1962 last_valid = value;
1963 }
1964 last_valid
1965 })
1966 .collect();
1967 stacked = Box::new(NormalizedTransmissionModel::new(
1968 stacked,
1969 &detector_energies,
1970 bi.anorm,
1971 bi.back_a,
1972 bi.back_b,
1973 bi.back_c,
1974 ));
1975 }
1976 } else {
1977 if let Some(bi) = bg_indices {
1981 stacked = Box::new(NormalizedTransmissionModel::new(
1982 stacked,
1983 config.energies(),
1984 bi.anorm,
1985 bi.back_a,
1986 bi.back_b,
1987 bi.back_c,
1988 ));
1989 }
1990 if let Some(bli) = bl_indices {
1991 stacked = Box::new(
1992 MultiplicativeBaselineModel::new(
1993 stacked,
1994 config.energies(),
1995 nereids_fitting::transmission_model::baseline_reference_energy_active(
1996 config.energies(),
1997 active_mask.as_deref(),
1998 ),
1999 bli.b0,
2000 bli.b1,
2001 bli.b2,
2002 )
2003 .with_active_mask(active_mask.as_deref()),
2005 );
2006 }
2007 }
2008 let background =
2018 (!detector_background.iter().all(|&v| v == 0.0)).then_some(detector_background);
2019 let objective = JointPoissonObjective {
2020 model: &*stacked,
2021 o: flux,
2022 s: sample_counts,
2023 c,
2024 active_mask: active_mask.as_deref(),
2025 open_background: background,
2026 sample_background: background,
2027 };
2028 let mut cfg = jp_cfg.clone();
2029 cfg.compute_covariance = config.compute_covariance;
2030 cfg.scale_by_chi2 = config.scale_by_chi2;
2031 let result = joint_poisson::joint_poisson_fit(&objective, &mut params, &cfg)
2039 .map_err(PipelineError::Fitting)?;
2040
2041 let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();
2043
2044 let (uncertainties, temperature_k_unc) = if let Some(ref unc_all) = result.uncertainties {
2045 let free_idx = params.free_indices();
2049 let dens_unc: Vec<f64> = (0..n_density_params)
2050 .map(|i| free_uncertainty(&free_idx, unc_all, i).unwrap_or(f64::NAN))
2051 .collect();
2052 let t_unc = temperature_index.and_then(|idx| free_uncertainty(&free_idx, unc_all, idx));
2053 (Some(dens_unc), t_unc)
2054 } else {
2055 (None, None)
2056 };
2057 let fitted_temp = temperature_index.map(|idx| result.params[idx]);
2058
2059 let converged = result.gn_converged || result.polish_converged;
2064
2065 let (anorm_out, bg_abc_out) = if let Some(bi) = bg_indices {
2070 (
2071 result.params[bi.anorm],
2072 [
2073 result.params[bi.back_a],
2074 result.params[bi.back_b],
2075 result.params[bi.back_c],
2076 ],
2077 )
2078 } else {
2079 (1.0, [0.0, 0.0, 0.0])
2080 };
2081
2082 let (baseline_out, baseline_e_ref_out) = if let Some(bli) = bl_indices {
2084 (
2085 Some([
2086 result.params[bli.b0],
2087 result.params[bli.b1],
2088 result.params[bli.b2],
2089 ]),
2090 Some(
2091 nereids_fitting::transmission_model::baseline_reference_energy_active(
2092 config.energies(),
2093 active_mask.as_deref(),
2094 ),
2095 ),
2096 )
2097 } else {
2098 (None, None)
2099 };
2100
2101 Ok(SpectrumFitResult {
2102 densities,
2103 uncertainties,
2104 reduced_chi_squared: result.deviance_per_dof,
2108 converged,
2109 iterations: result.gn_iterations + result.polish_iterations,
2110 temperature_k: fitted_temp,
2111 temperature_k_unc,
2112 anorm: anorm_out,
2113 background: bg_abc_out,
2114 back_d: None,
2119 back_f: None,
2120 t0_us: energy_scale_indices.map(|(t0_idx, _)| result.params[t0_idx]),
2121 l_scale: energy_scale_indices.map(|(_, ls_idx)| result.params[ls_idx]),
2122 energy_scale_flight_path_m: energy_scale_indices.map(|_| config.flight_path_m),
2123 deviance_per_dof: Some(result.deviance_per_dof),
2124 baseline: baseline_out,
2125 baseline_e_ref_ev: baseline_e_ref_out,
2126 warnings: degenerate_normalization_warning(config)
2127 .into_iter()
2128 .collect(),
2129 })
2130}
2131
2132fn build_density_params(config: &UnifiedFitConfig) -> Vec<FitParameter> {
2135 config
2136 .initial_densities
2137 .iter()
2138 .enumerate()
2139 .map(|(i, &d)| {
2140 let name = config
2141 .isotope_names
2142 .get(i)
2143 .cloned()
2144 .unwrap_or_else(|| format!("isotope_{i}"));
2145 if config.density_is_fixed(i) {
2150 FitParameter::fixed(name, d)
2151 } else {
2152 FitParameter::non_negative(name, d)
2153 }
2154 })
2155 .collect()
2156}
2157
2158pub(crate) const TEMPERATURE_BOUNDS_K: (f64, f64) = (1.0, 5000.0);
2166
2167fn append_temperature_param(
2172 param_vec: &mut Vec<FitParameter>,
2173 config: &UnifiedFitConfig,
2174) -> Option<usize> {
2175 if !config.fit_temperature {
2176 return None;
2177 }
2178 let idx = param_vec.len();
2179 param_vec.push(FitParameter {
2180 name: "temperature_k".into(),
2181 value: config.temperature_k,
2182 lower: TEMPERATURE_BOUNDS_K.0,
2183 upper: TEMPERATURE_BOUNDS_K.1,
2184 fixed: false,
2185 });
2186 Some(idx)
2187}
2188
2189fn detect_transmission_dips(measured: &[f64], energies: &[f64]) -> Vec<(f64, f64)> {
2195 let n = measured.len();
2196 if n < 5 || energies.len() != n {
2197 return Vec::new();
2198 }
2199 let mut sm = measured.to_vec();
2200 for i in 1..n - 1 {
2201 sm[i] = (measured[i - 1] + measured[i] + measured[i + 1]) / 3.0;
2202 }
2203 let mut sorted = sm.clone();
2205 sorted.sort_by(f64::total_cmp);
2206 let baseline = sorted[((0.9 * (n as f64 - 1.0)).round() as usize).min(n - 1)];
2207 let mut dips: Vec<(f64, f64)> = Vec::new();
2208 let mut max_depth = 0.0f64;
2209 for i in 1..n - 1 {
2210 if sm[i] < sm[i - 1] && sm[i] <= sm[i + 1] {
2211 let depth = baseline - sm[i];
2212 if depth > 0.0 {
2213 dips.push((energies[i], depth));
2214 max_depth = max_depth.max(depth);
2215 }
2216 }
2217 }
2218 dips.retain(|&(_, d)| d >= 0.2 * max_depth);
2221 dips
2222}
2223
2224fn peak_match_energy_scale_seed(
2240 measured: &[f64],
2241 energies: &[f64],
2242 config: &UnifiedFitConfig,
2243 flight_path_m: f64,
2244 t0_bounds: (f64, f64),
2245 l_scale_bounds: (f64, f64),
2246) -> Option<(f64, f64)> {
2247 let n = energies.len();
2248 if n < 5 || measured.len() != n {
2249 return None;
2250 }
2251 let refs: Vec<&ResonanceData> = config.resonance_data().iter().collect();
2253 let (e_lo, e_hi) = (
2254 energies[0].min(energies[n - 1]),
2255 energies[0].max(energies[n - 1]),
2256 );
2257 let res_e: Vec<f64> = nereids_physics::transmission::resonance_center_energies(&refs)
2258 .into_iter()
2259 .filter(|&e| e > e_lo && e < e_hi)
2260 .collect();
2261 if res_e.len() < 2 {
2262 return None;
2263 }
2264 let dips = detect_transmission_dips(measured, energies);
2265 if dips.len() < 2 {
2266 return None;
2267 }
2268 let min_spacing = res_e
2283 .windows(2)
2284 .map(|w| (w[1] - w[0]).abs())
2285 .fold(f64::INFINITY, f64::min);
2286 let mut grid_steps: Vec<f64> = energies.windows(2).map(|w| (w[1] - w[0]).abs()).collect();
2287 grid_steps.sort_by(f64::total_cmp);
2288 let grid_res = grid_steps.get(grid_steps.len() / 2).copied().unwrap_or(0.0);
2289 let match_tol = (0.5 * min_spacing).max(grid_res);
2290 let mut best: Vec<Option<(f64, f64)>> = vec![None; res_e.len()];
2291 for &(e_dip, depth) in &dips {
2292 let Some((k, &re)) = res_e
2293 .iter()
2294 .enumerate()
2295 .min_by(|(_, a), (_, b)| (*a - e_dip).abs().total_cmp(&(*b - e_dip).abs()))
2296 else {
2297 continue;
2298 };
2299 if (re - e_dip).abs() > match_tol {
2300 continue;
2301 }
2302 match best[k] {
2303 Some((_, d)) if d >= depth => {}
2304 _ => best[k] = Some((e_dip, depth)),
2305 }
2306 }
2307 let tof_factor = (0.5 * NEUTRON_MASS_KG / EV_TO_JOULES).sqrt() * 1.0e6;
2311 let c = tof_factor * flight_path_m;
2312 let pairs: Vec<(f64, f64)> = res_e
2313 .iter()
2314 .zip(best.iter())
2315 .filter_map(|(&re, b)| b.map(|(e_dip, _)| (c / re.sqrt(), c / e_dip.sqrt())))
2316 .collect();
2317 if pairs.len() < 2 {
2318 return None;
2319 }
2320 let m = pairs.len() as f64;
2322 let mean_x = pairs.iter().map(|p| p.0).sum::<f64>() / m;
2323 let mean_y = pairs.iter().map(|p| p.1).sum::<f64>() / m;
2324 let mut sxx = 0.0f64;
2325 let mut sxy = 0.0f64;
2326 for &(x, y) in &pairs {
2327 sxx += (x - mean_x) * (x - mean_x);
2328 sxy += (x - mean_x) * (y - mean_y);
2329 }
2330 if sxx <= 0.0 {
2331 return None;
2332 }
2333 let l_scale = sxy / sxx;
2334 let t0 = mean_y - l_scale * mean_x;
2335 if !t0.is_finite() || !l_scale.is_finite() {
2336 return None;
2337 }
2338 if t0 < t0_bounds.0
2343 || t0 > t0_bounds.1
2344 || l_scale < l_scale_bounds.0
2345 || l_scale > l_scale_bounds.1
2346 {
2347 return None;
2348 }
2349 Some((t0, l_scale))
2350}
2351
2352fn seed_energy_scale_in_params(
2356 param_vec: &mut [FitParameter],
2357 energy_scale_indices: Option<(usize, usize)>,
2358 measured_transmission: &[f64],
2359 config: &UnifiedFitConfig,
2360) {
2361 if !config.energy_scale_seed_enabled {
2364 return;
2365 }
2366 let Some((t0_idx, ls_idx)) = energy_scale_indices else {
2367 return;
2368 };
2369 let t0_b = (param_vec[t0_idx].lower, param_vec[t0_idx].upper);
2370 let ls_b = (param_vec[ls_idx].lower, param_vec[ls_idx].upper);
2371 if let Some((t0_seed, ls_seed)) = peak_match_energy_scale_seed(
2372 measured_transmission,
2373 config.energies(),
2374 config,
2375 config.flight_path_m,
2376 t0_b,
2377 ls_b,
2378 ) {
2379 param_vec[t0_idx].value = t0_seed;
2380 param_vec[ls_idx].value = ls_seed;
2381 }
2382}
2383
2384const ENERGY_SCALE_T0_BOUND_US: f64 = 10.0;
2388const ENERGY_SCALE_L_SCALE_LO: f64 = 0.99;
2392const ENERGY_SCALE_L_SCALE_HI: f64 = 1.01;
2394
2395fn append_energy_scale_params(
2401 param_vec: &mut Vec<FitParameter>,
2402 config: &UnifiedFitConfig,
2403) -> Option<(usize, usize)> {
2404 if !config.fit_energy_scale {
2405 return None;
2406 }
2407 let t0_idx = param_vec.len();
2408 param_vec.push(FitParameter {
2409 name: "t0_us".into(),
2410 value: config.t0_init_us,
2411 lower: -ENERGY_SCALE_T0_BOUND_US,
2412 upper: ENERGY_SCALE_T0_BOUND_US,
2413 fixed: false,
2414 });
2415 let ls_idx = param_vec.len();
2416 param_vec.push(FitParameter {
2417 name: "l_scale".into(),
2418 value: config.l_scale_init,
2419 lower: ENERGY_SCALE_L_SCALE_LO,
2420 upper: ENERGY_SCALE_L_SCALE_HI,
2421 fixed: false,
2422 });
2423 Some((t0_idx, ls_idx))
2424}
2425
2426pub(crate) fn validate_transmission_background(bg: &BackgroundConfig) -> Result<(), PipelineError> {
2441 if bg.fit_back_d != bg.fit_back_f {
2442 return Err(PipelineError::InvalidParameter(format!(
2443 "transmission_background: fit_back_d ({}) and fit_back_f ({}) \
2444 must both be true or both be false. The exponential tail \
2445 wrapper (BackD · exp(−BackF / √E)) requires both parameters \
2446 together; enabling only one leaves the other registered but \
2447 unused, silently producing the initial value as the fitted \
2448 result. Either enable both (to fit the exponential tail) or \
2449 disable both (4-term wrapper without exponential).",
2450 bg.fit_back_d, bg.fit_back_f,
2451 )));
2452 }
2453 Ok(())
2454}
2455
2456pub(crate) fn validate_multiplicative_baseline(
2472 config: &UnifiedFitConfig,
2473) -> Result<(), PipelineError> {
2474 let Some(bl) = config.multiplicative_baseline() else {
2475 return Ok(());
2476 };
2477 for (name, init, (lo, hi)) in [
2478 ("b0", bl.b0_init, bl.b0_bounds),
2479 ("b1", bl.b1_init, bl.b1_bounds),
2480 ("b2", bl.b2_init, bl.b2_bounds),
2481 ] {
2482 if !init.is_finite() {
2483 return Err(PipelineError::InvalidParameter(format!(
2484 "multiplicative baseline: {name}_init must be finite, got {init}"
2485 )));
2486 }
2487 if !lo.is_finite() || !hi.is_finite() || lo >= hi {
2488 return Err(PipelineError::InvalidParameter(format!(
2489 "multiplicative baseline: {name}_bounds must be finite with \
2490 lower < upper, got ({lo}, {hi})"
2491 )));
2492 }
2493 if init < lo || init > hi {
2494 return Err(PipelineError::InvalidParameter(format!(
2495 "multiplicative baseline: {name}_init = {init} lies outside \
2496 {name}_bounds ({lo}, {hi})"
2497 )));
2498 }
2499 }
2500 let e_ref = config.baseline_reference_energy();
2513 if !e_ref.is_finite() || e_ref <= 0.0 {
2514 return Err(PipelineError::InvalidParameter(format!(
2515 "multiplicative baseline: reference energy sqrt(E_min*E_max) is \
2516 invalid ({e_ref}) — check the energy grid"
2517 )));
2518 }
2519 let active_mask = nereids_fitting::active_mask::build_active_mask(
2520 config.energies(),
2521 config.fit_energy_range(),
2522 );
2523 for (i, &e) in config.energies().iter().enumerate() {
2524 if active_mask.as_ref().is_some_and(|m| !m[i]) {
2525 continue;
2526 }
2527 let z = (e / e_ref).ln();
2528 let b = bl.b0_init + bl.b1_init * z + bl.b2_init * z * z;
2529 let positive = b.is_finite() && b > 0.0;
2530 if !positive {
2531 return Err(PipelineError::InvalidParameter(format!(
2532 "multiplicative baseline: initial B(E) = {b} is not strictly \
2533 positive at E = {e} eV (inside the fit window) — adjust the \
2534 b0/b1/b2 inits"
2535 )));
2536 }
2537 }
2538 if config
2539 .transmission_background()
2540 .is_some_and(|bg| bg.fit_anorm)
2541 {
2542 return Err(PipelineError::InvalidParameter(
2543 "multiplicative baseline and a FREE SAMMY Anorm cannot be fitted \
2544 together: b0 and Anorm are degenerate normalizations. Set \
2545 BackgroundConfig::fit_anorm = false to combine the additive ABC \
2546 background with the baseline."
2547 .into(),
2548 ));
2549 }
2550 Ok(())
2551}
2552
2553pub(crate) fn degenerate_normalization_warning(config: &UnifiedFitConfig) -> Option<String> {
2559 let anorm_free = config
2560 .transmission_background()
2561 .is_some_and(|bg| bg.fit_anorm);
2562 if anorm_free && config.fit_temperature && config.n_free_density_params() >= 1 {
2563 Some(
2564 "fit configuration frees Anorm, temperature, AND at least one \
2565 density together — a degenerate normalization trio on real data \
2566 (observed: T ran to 4471 K with chi2/nu 932 and no warning). \
2567 Consider with_multiplicative_baseline (bounded normalization) \
2568 and/or with_fix_densities (known areal density)."
2569 .to_string(),
2570 )
2571 } else {
2572 None
2573 }
2574}
2575
2576pub(crate) fn validate_precomputed_cross_sections(
2577 config: &UnifiedFitConfig,
2578) -> Result<(), PipelineError> {
2579 let Some(xs) = config.precomputed_cross_sections() else {
2580 return Ok(());
2581 };
2582 if xs.sigma.is_empty() {
2583 return Err(PipelineError::ShapeMismatch(
2584 "precomputed_cross_sections must not be empty".into(),
2585 ));
2586 }
2587 let n_work = xs.layout.energies.len();
2588 for (i, row) in xs.sigma.iter().enumerate() {
2589 if row.len() != n_work {
2590 return Err(PipelineError::ShapeMismatch(format!(
2591 "precomputed_cross_sections row {i} has length {} but its grid has {n_work} \
2592 energies",
2593 row.len(),
2594 )));
2595 }
2596 if let Some(j) = row.iter().position(|s| !s.is_finite()) {
2597 return Err(PipelineError::ShapeMismatch(format!(
2598 "precomputed_cross_sections row {i} has non-finite σ at energy index {j}: {}",
2599 row[j],
2600 )));
2601 }
2602 }
2603
2604 let n_params = config.n_density_params();
2605 let member_rows = config.density_indices.as_ref().map(|di| di.len());
2606 let row_ok = xs.sigma.len() == n_params || member_rows == Some(xs.sigma.len());
2607 if !row_ok {
2608 let expected = match member_rows {
2609 Some(m) if m != n_params => format!("{n_params} (collapsed) or {m} (per-member)"),
2610 _ => format!("{n_params}"),
2611 };
2612 return Err(PipelineError::ShapeMismatch(format!(
2613 "precomputed_cross_sections has {} rows but expected {expected}",
2614 xs.sigma.len(),
2615 )));
2616 }
2617
2618 let instrument = if config.exact_count_response().is_some() {
2619 None
2620 } else {
2621 config.resolution().map(|r| InstrumentParams {
2622 resolution: r.clone(),
2623 })
2624 };
2625 let rd_refs: Vec<&ResonanceData> = config.resonance_data().iter().collect();
2626 let expected = nereids_physics::transmission::resolution_working_grid(
2627 config.energies(),
2628 instrument.as_ref(),
2629 &rd_refs,
2630 )
2631 .map_err(PipelineError::Transmission)?;
2632 let same_energies = xs.layout.energies.len() == expected.energies.len()
2633 && xs
2634 .layout
2635 .energies
2636 .iter()
2637 .zip(&expected.energies)
2638 .all(|(a, b)| a.to_bits() == b.to_bits());
2639 if !same_energies || xs.layout.data_indices != expected.data_indices {
2640 return Err(PipelineError::ShapeMismatch(format!(
2641 "precomputed_cross_sections layout ({} energies, {} data points) is not the \
2642 working grid this fit broadens on ({} energies, {} data points)",
2643 xs.layout.energies.len(),
2644 xs.layout.data_indices.len(),
2645 expected.energies.len(),
2646 expected.data_indices.len(),
2647 )));
2648 }
2649 Ok(())
2650}
2651
2652pub(crate) fn count_free_params(config: &UnifiedFitConfig) -> usize {
2682 let mut n_free = config.n_free_density_params();
2686 if config.fit_temperature {
2687 n_free += 1;
2688 }
2689 if config.fit_energy_scale {
2690 n_free += 2;
2691 }
2692 if let Some(bg) = config.transmission_background.as_ref() {
2693 n_free += usize::from(bg.fit_anorm);
2694 n_free += usize::from(bg.fit_back_a);
2695 n_free += usize::from(bg.fit_back_b);
2696 n_free += usize::from(bg.fit_back_c);
2697 n_free += usize::from(bg.fit_back_d);
2698 n_free += usize::from(bg.fit_back_f);
2699 }
2700 if let Some(bl) = config.multiplicative_baseline.as_ref() {
2701 n_free += usize::from(bl.fit_b0);
2702 n_free += usize::from(bl.fit_b1);
2703 n_free += usize::from(bl.fit_b2);
2704 }
2705 n_free
2706}
2707
2708pub(crate) fn required_active_bins(config: &UnifiedFitConfig) -> usize {
2727 count_free_params(config).max(2)
2728}
2729
2730fn append_background_params(
2731 param_vec: &mut Vec<FitParameter>,
2732 bg: &BackgroundConfig,
2733) -> BackgroundIndices {
2734 let anorm = param_vec.len();
2739 param_vec.push(if bg.fit_anorm {
2740 FitParameter {
2741 name: "anorm".into(),
2742 value: bg.anorm_init,
2743 lower: 0.5,
2744 upper: 2.0,
2745 fixed: false,
2746 }
2747 } else {
2748 FitParameter::fixed("anorm", bg.anorm_init)
2749 });
2750 let back_a = param_vec.len();
2756 param_vec.push(if bg.fit_back_a {
2757 FitParameter {
2758 name: "back_a".into(),
2759 value: bg.back_a_init,
2760 lower: -0.5,
2761 upper: 0.5,
2762 fixed: false,
2763 }
2764 } else {
2765 FitParameter::fixed("back_a", bg.back_a_init)
2766 });
2767 let back_b = param_vec.len();
2768 param_vec.push(if bg.fit_back_b {
2769 FitParameter {
2770 name: "back_b".into(),
2771 value: bg.back_b_init,
2772 lower: -0.5,
2773 upper: 0.5,
2774 fixed: false,
2775 }
2776 } else {
2777 FitParameter::fixed("back_b", bg.back_b_init)
2778 });
2779 let back_c = param_vec.len();
2780 param_vec.push(if bg.fit_back_c {
2781 FitParameter {
2782 name: "back_c".into(),
2783 value: bg.back_c_init,
2784 lower: -0.5,
2785 upper: 0.5,
2786 fixed: false,
2787 }
2788 } else {
2789 FitParameter::fixed("back_c", bg.back_c_init)
2790 });
2791
2792 let back_d = if bg.fit_back_d {
2800 let idx = param_vec.len();
2801 param_vec.push(FitParameter {
2802 name: "back_d".into(),
2803 value: bg.back_d_init,
2804 lower: 0.0,
2805 upper: 1.0,
2806 fixed: false,
2807 });
2808 Some(idx)
2809 } else {
2810 None
2811 };
2812 let back_f = if bg.fit_back_f {
2813 let idx = param_vec.len();
2814 param_vec.push(FitParameter {
2815 name: "back_f".into(),
2816 value: bg.back_f_init,
2817 lower: 0.0,
2818 upper: 100.0,
2819 fixed: false,
2820 });
2821 Some(idx)
2822 } else {
2823 None
2824 };
2825
2826 BackgroundIndices {
2827 anorm,
2828 back_a,
2829 back_b,
2830 back_c,
2831 back_d,
2832 back_f,
2833 }
2834}
2835
2836fn append_multiplicative_baseline_params(
2848 param_vec: &mut Vec<FitParameter>,
2849 bl: &MultiplicativeBaselineConfig,
2850) -> BaselineIndices {
2851 let b0 = param_vec.len();
2852 param_vec.push(if bl.fit_b0 {
2853 FitParameter {
2854 name: "baseline_b0".into(),
2855 value: bl.b0_init,
2856 lower: bl.b0_bounds.0,
2857 upper: bl.b0_bounds.1,
2858 fixed: false,
2859 }
2860 } else {
2861 FitParameter::fixed("baseline_b0", bl.b0_init)
2862 });
2863 let b1 = param_vec.len();
2864 param_vec.push(if bl.fit_b1 {
2865 FitParameter {
2866 name: "baseline_b1".into(),
2867 value: bl.b1_init,
2868 lower: bl.b1_bounds.0,
2869 upper: bl.b1_bounds.1,
2870 fixed: false,
2871 }
2872 } else {
2873 FitParameter::fixed("baseline_b1", bl.b1_init)
2874 });
2875 let b2 = param_vec.len();
2876 param_vec.push(if bl.fit_b2 {
2877 FitParameter {
2878 name: "baseline_b2".into(),
2879 value: bl.b2_init,
2880 lower: bl.b2_bounds.0,
2881 upper: bl.b2_bounds.1,
2882 fixed: false,
2883 }
2884 } else {
2885 FitParameter::fixed("baseline_b2", bl.b2_init)
2886 });
2887
2888 BaselineIndices { b0, b1, b2 }
2889}
2890
2891fn build_energy_scale_transmission_model(
2914 config: &UnifiedFitConfig,
2915 t0_idx: usize,
2916 ls_idx: usize,
2917 temperature_index: Option<usize>,
2918) -> Result<Box<dyn FitModel>, PipelineError> {
2919 let instrument = config
2920 .resolution
2921 .clone()
2922 .map(|r| Arc::new(InstrumentParams { resolution: r }));
2923 let n_iso = config.resonance_data.len();
2929 let density_indices = config
2930 .density_indices
2931 .clone()
2932 .unwrap_or_else(|| (0..n_iso).collect());
2933 let density_ratios = config
2934 .density_ratios
2935 .clone()
2936 .unwrap_or_else(|| vec![1.0; n_iso]);
2937 let mut es_model = EnergyScaleTransmissionModel::new(
2938 Arc::new(config.resonance_data.clone()),
2939 Arc::new(density_indices),
2940 Arc::new(density_ratios),
2941 config.temperature_k,
2942 config.energies.clone(),
2943 config.flight_path_m,
2944 t0_idx,
2945 ls_idx,
2946 instrument,
2947 )
2948 .with_temperature_index(temperature_index)
2952 .map_err(|e| PipelineError::InvalidParameter(format!("energy-scale model: {e}")))?;
2953 if let Some(method) = config.tzero_jacobian_method {
2954 es_model = es_model.with_jacobian_method(method);
2955 }
2956 Ok(Box::new(es_model))
2957}
2958
2959fn build_transmission_model(
2961 config: &UnifiedFitConfig,
2962 n_density_params: usize,
2963 temperature_index: Option<usize>,
2964) -> Result<Box<dyn FitModel>, PipelineError> {
2965 let n_params = config.n_density_params();
2966 let instrument = config
2967 .resolution
2968 .clone()
2969 .map(|r| Arc::new(InstrumentParams { resolution: r }));
2970
2971 if !config.fit_temperature {
2972 let xs = match &config.precomputed_cross_sections {
2973 Some(xs) => xs.clone(),
2974 None => PrecomputedXs::from(
2975 nereids_physics::transmission::broadened_cross_sections_on_working_grid(
2976 config.energies(),
2977 &config.resonance_data,
2978 config.temperature_k,
2979 instrument.as_deref(),
2980 None,
2981 )
2982 .map_err(PipelineError::Transmission)?,
2983 ),
2984 };
2985 let cross_sections = match (&config.density_indices, &config.density_ratios) {
2986 (Some(di), Some(dr)) if xs.sigma.len() == di.len() && di.len() == dr.len() => {
2987 let mut eff = vec![vec![0.0f64; xs.sigma[0].len()]; n_params];
2988 for ((&idx, &ratio), member) in di.iter().zip(dr.iter()).zip(xs.sigma.iter()) {
2989 for (j, &sigma) in member.iter().enumerate() {
2990 eff[idx][j] += ratio * sigma;
2991 }
2992 }
2993 Arc::new(eff)
2994 }
2995 _ => Arc::clone(&xs.sigma),
2996 };
2997 let (resolution_plan, sparse_cubature_plan, sparse_scalar_plan) = if instrument.is_some() {
2998 (
2999 config.precomputed_resolution_plan.clone(),
3000 config.precomputed_sparse_cubature_plan.clone(),
3001 config.precomputed_sparse_scalar_plan.clone(),
3002 )
3003 } else {
3004 (None, None, None)
3005 };
3006 return Ok(Box::new(PrecomputedTransmissionModel {
3007 cross_sections,
3008 density_indices: Arc::new((0..n_params).collect()),
3009 instrument,
3010 resolution_plan,
3011 sparse_cubature_plan,
3012 sparse_scalar_plan,
3013 layout: Arc::clone(&xs.layout),
3014 }));
3015 }
3016
3017 let base_xs = config.precomputed_base_xs.clone();
3018 let density_ratios = config
3019 .density_ratios
3020 .clone()
3021 .unwrap_or_else(|| vec![1.0; n_density_params]);
3022 let density_indices = config
3023 .density_indices
3024 .clone()
3025 .unwrap_or_else(|| (0..n_density_params).collect());
3026 let resolution_plan = if instrument.is_some() {
3027 config.precomputed_resolution_plan.clone()
3028 } else {
3029 None
3030 };
3031 let sparse_cubature_plan = if instrument.is_some() {
3032 config.precomputed_sparse_cubature_plan.clone()
3033 } else {
3034 None
3035 };
3036 let sparse_scalar_plan = if instrument.is_some() {
3037 config.precomputed_sparse_scalar_plan.clone()
3038 } else {
3039 None
3040 };
3041 Ok(Box::new(
3042 TransmissionFitModel::new(
3043 config.energies.clone(),
3044 config.resonance_data.clone(),
3045 config.temperature_k,
3046 instrument,
3047 (density_indices, density_ratios),
3048 temperature_index,
3049 base_xs,
3050 )?
3051 .with_resolution_plan(resolution_plan)
3052 .with_sparse_cubature_plan(sparse_cubature_plan)
3053 .with_sparse_scalar_plan(sparse_scalar_plan),
3054 ))
3055}
3056
3057fn free_uncertainty(free_indices: &[usize], unc_all: &[f64], full_index: usize) -> Option<f64> {
3069 free_indices
3070 .iter()
3071 .position(|&fi| fi == full_index)
3072 .and_then(|pos| unc_all.get(pos).copied())
3073}
3074
3075fn extract_result(
3077 config: &UnifiedFitConfig,
3078 result: &LmResult,
3079 n_density_params: usize,
3080 free_indices: &[usize],
3081 bg_indices: Option<BackgroundIndices>,
3082 bl_indices: Option<BaselineIndices>,
3083) -> Result<SpectrumFitResult, PipelineError> {
3084 let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();
3085
3086 let (baseline, baseline_e_ref_ev) = if let Some(bli) = bl_indices {
3090 (
3091 Some([
3092 result.params[bli.b0],
3093 result.params[bli.b1],
3094 result.params[bli.b2],
3095 ]),
3096 Some(config.baseline_reference_energy()),
3097 )
3098 } else {
3099 (None, None)
3100 };
3101
3102 let (anorm, background, back_d, back_f): (f64, [f64; 3], Option<f64>, Option<f64>) =
3103 if let Some(bi) = bg_indices {
3104 let bd = bi.back_d.map(|i| result.params[i]);
3107 let bf = bi.back_f.map(|i| result.params[i]);
3108 (
3109 result.params[bi.anorm],
3110 [
3111 result.params[bi.back_a],
3112 result.params[bi.back_b],
3113 result.params[bi.back_c],
3114 ],
3115 bd,
3116 bf,
3117 )
3118 } else {
3119 (1.0, [0.0, 0.0, 0.0], None, None)
3120 };
3121
3122 let (uncertainties, temperature_k, temperature_k_unc) = if result.converged {
3123 match &result.uncertainties {
3124 Some(unc_all) => {
3125 let (temp_k, temp_unc) = if config.fit_temperature {
3134 (
3135 Some(result.params[n_density_params]),
3136 Some(
3137 free_uncertainty(free_indices, unc_all, n_density_params)
3138 .unwrap_or(f64::NAN),
3139 ),
3140 )
3141 } else {
3142 (None, None)
3143 };
3144 let unc: Vec<f64> = (0..n_density_params)
3145 .map(|i| free_uncertainty(free_indices, unc_all, i).unwrap_or(f64::NAN))
3146 .collect();
3147 (Some(unc), temp_k, temp_unc)
3148 }
3149 None => {
3150 let temp_k = if config.fit_temperature {
3151 Some(result.params[n_density_params])
3152 } else {
3153 None
3154 };
3155 (None, temp_k, None)
3156 }
3157 }
3158 } else {
3159 let temp_k = if config.fit_temperature {
3160 Some(result.params[n_density_params])
3161 } else {
3162 None
3163 };
3164 (None, temp_k, None)
3165 };
3166
3167 Ok(SpectrumFitResult {
3168 densities,
3169 uncertainties,
3170 reduced_chi_squared: result.reduced_chi_squared,
3171 converged: result.converged,
3172 iterations: result.iterations,
3173 temperature_k,
3174 temperature_k_unc,
3175 anorm,
3176 background,
3177 back_d,
3178 back_f,
3179 t0_us: None,
3180 l_scale: None,
3181 energy_scale_flight_path_m: None,
3182 deviance_per_dof: None,
3183 baseline,
3184 baseline_e_ref_ev,
3185 warnings: degenerate_normalization_warning(config)
3186 .into_iter()
3187 .collect(),
3188 })
3189}
3190
3191pub struct ModelJacobianResult {
3199 pub jacobian: lm::FlatMatrix,
3201 pub fisher: lm::FlatMatrix,
3203 pub model_prediction: Vec<f64>,
3205 pub param_names: Vec<String>,
3207}
3208
3209pub fn evaluate_jacobian_and_fisher(
3210 config: &UnifiedFitConfig,
3211 flux: &[f64],
3212 background: &[f64],
3213) -> Result<ModelJacobianResult, PipelineError> {
3214 if config.exact_count_response().is_some() || config.resolution().is_some() {
3222 return Err(PipelineError::InvalidParameter(
3223 "evaluate_jacobian_and_fisher does not implement the exact two-arm \
3224 detector operator required for resolved counts (the separate-arm \
3225 model R[Phi] and R[Phi*T]); use fit_counts_spectrum_typed with \
3226 exact_count_response for resolved counts, or drop the instrument \
3227 resolution for this research helper"
3228 .into(),
3229 ));
3230 }
3231 validate_counts_resolution_route(true, flux.len(), config)?;
3234
3235 validate_precomputed_cross_sections(config)?;
3244
3245 if config.multiplicative_baseline.is_some() {
3253 return Err(PipelineError::InvalidParameter(
3254 "evaluate_jacobian_and_fisher does not support a multiplicative \
3255 baseline: its parameter layout is frozen for research-script \
3256 compatibility (kl_b0/kl_b1 background stand-ins). Remove \
3257 with_multiplicative_baseline from the config for this helper."
3258 .into(),
3259 ));
3260 }
3261
3262 let n_density_params = config.n_density_params();
3263
3264 let mut param_vec = build_density_params(config);
3271
3272 let temperature_index = if config.fit_temperature {
3273 let idx = param_vec.len();
3274 param_vec.push(FitParameter {
3275 name: "temperature_k".into(),
3276 value: config.temperature_k,
3277 lower: TEMPERATURE_BOUNDS_K.0,
3278 upper: TEMPERATURE_BOUNDS_K.1,
3279 fixed: false,
3280 });
3281 Some(idx)
3282 } else {
3283 None
3284 };
3285
3286 let kl_bg = if config.transmission_background.is_some() {
3287 let base = param_vec.len();
3288 param_vec.push(FitParameter {
3289 name: "kl_b0".into(),
3290 value: 0.0,
3291 lower: 0.0,
3292 upper: 0.5,
3293 fixed: false,
3294 });
3295 param_vec.push(FitParameter {
3296 name: "kl_b1".into(),
3297 value: 0.0,
3298 lower: 0.0,
3299 upper: 0.5,
3300 fixed: false,
3301 });
3302 Some((base, base + 1))
3303 } else {
3304 None
3305 };
3306
3307 let counts_bg = if let Some(bg) = config.counts_background() {
3308 let alpha1_idx = param_vec.len();
3309 param_vec.push(if bg.fit_alpha_1 {
3310 FitParameter {
3311 name: "alpha_1".into(),
3312 value: bg.alpha_1_init,
3313 lower: 0.0,
3314 upper: 10.0,
3315 fixed: false,
3316 }
3317 } else {
3318 FitParameter::fixed("alpha_1", bg.alpha_1_init)
3319 });
3320 let alpha2_idx = param_vec.len();
3321 param_vec.push(if bg.fit_alpha_2 {
3322 FitParameter {
3323 name: "alpha_2".into(),
3324 value: bg.alpha_2_init,
3325 lower: 0.0,
3326 upper: 10.0,
3327 fixed: false,
3328 }
3329 } else {
3330 FitParameter::fixed("alpha_2", bg.alpha_2_init)
3331 });
3332 Some((alpha1_idx, alpha2_idx))
3333 } else {
3334 None
3335 };
3336
3337 let params = ParameterSet::new(param_vec);
3338 let all_vals = params.all_values();
3339 let free_idx = params.free_indices();
3340 let n_free = free_idx.len();
3341
3342 let param_names: Vec<String> = free_idx
3344 .iter()
3345 .map(|&i| params.params[i].name.to_string())
3346 .collect();
3347
3348 let t_model = build_transmission_model(config, n_density_params, temperature_index)?;
3349
3350 let evaluate_and_jacobian =
3353 |model: &dyn FitModel| -> Result<(Vec<f64>, lm::FlatMatrix), PipelineError> {
3354 let y_model = model.evaluate(&all_vals)?;
3355 let jac = model
3356 .analytical_jacobian(&all_vals, &free_idx, &y_model)
3357 .ok_or_else(|| {
3358 PipelineError::InvalidParameter(
3359 "analytical Jacobian not available for this model configuration".into(),
3360 )
3361 })?;
3362 Ok((y_model, jac))
3363 };
3364
3365 let (y_model, jac) = if let Some((b0_idx, b1_idx)) = kl_bg {
3366 let inv_sqrt_e: Vec<f64> = config
3367 .energies()
3368 .iter()
3369 .map(|&e| 1.0 / e.max(1e-10).sqrt())
3370 .collect();
3371 let wrapped = poisson::TransmissionKLBackgroundModel {
3372 inner: &*t_model,
3373 inv_sqrt_energies: inv_sqrt_e,
3374 b0_index: b0_idx,
3375 b1_index: b1_idx,
3376 n_params: params.params.len(),
3377 };
3378 if let Some((a1, a2)) = counts_bg {
3379 let cm = poisson::CountsBackgroundScaleModel {
3380 transmission_model: &wrapped,
3381 flux,
3382 background,
3383 alpha1_index: a1,
3384 alpha2_index: a2,
3385 n_params: params.params.len(),
3386 };
3387 evaluate_and_jacobian(&cm)?
3388 } else {
3389 let cm = poisson::CountsModel {
3390 transmission_model: &wrapped,
3391 flux,
3392 background,
3393 n_params: params.params.len(),
3394 };
3395 evaluate_and_jacobian(&cm)?
3396 }
3397 } else if let Some((a1, a2)) = counts_bg {
3398 let cm = poisson::CountsBackgroundScaleModel {
3399 transmission_model: &*t_model,
3400 flux,
3401 background,
3402 alpha1_index: a1,
3403 alpha2_index: a2,
3404 n_params: params.params.len(),
3405 };
3406 evaluate_and_jacobian(&cm)?
3407 } else {
3408 let cm = poisson::CountsModel {
3409 transmission_model: &*t_model,
3410 flux,
3411 background,
3412 n_params: params.params.len(),
3413 };
3414 evaluate_and_jacobian(&cm)?
3415 };
3416
3417 let mut fisher = lm::FlatMatrix::zeros(n_free, n_free);
3419 for (i, &mu_i) in y_model.iter().enumerate() {
3420 let mu_inv = 1.0 / mu_i.max(1e-30);
3421 for a in 0..n_free {
3422 let ja = jac.get(i, a);
3423 for b in 0..=a {
3424 let jb = jac.get(i, b);
3425 *fisher.get_mut(a, b) += ja * jb * mu_inv;
3426 if a != b {
3427 *fisher.get_mut(b, a) += ja * jb * mu_inv;
3428 }
3429 }
3430 }
3431 }
3432
3433 Ok(ModelJacobianResult {
3434 jacobian: jac,
3435 fisher,
3436 model_prediction: y_model,
3437 param_names,
3438 })
3439}
3440
3441#[derive(Debug, PartialEq)]
3445pub enum FitConfigError {
3446 EmptyEnergies,
3448 EmptyResonanceData,
3450 DensityCountMismatch { densities: usize, isotopes: usize },
3452 NameCountMismatch { names: usize, isotopes: usize },
3454 GroupMemberCountMismatch {
3456 group_name: String,
3457 rd_count: usize,
3458 member_count: usize,
3459 },
3460 GroupMemberIsotopeMismatch {
3462 group_name: String,
3463 expected_z: u32,
3464 expected_a: u32,
3465 got_z: u32,
3466 got_a: u32,
3467 },
3468 NonFiniteTemperature(f64),
3470 NegativeTemperature(f64),
3472 InvalidFitEnergyRange(&'static str),
3474 DensityFreezeBeforeGroups,
3479}
3480
3481impl fmt::Display for FitConfigError {
3482 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
3483 match self {
3484 Self::EmptyEnergies => write!(f, "energy grid must be non-empty"),
3485 Self::EmptyResonanceData => write!(f, "resonance_data must be non-empty"),
3486 Self::DensityCountMismatch {
3487 densities,
3488 isotopes,
3489 } => write!(
3490 f,
3491 "initial_densities length ({densities}) must match number of density parameters ({isotopes})"
3492 ),
3493 Self::NameCountMismatch { names, isotopes } => write!(
3494 f,
3495 "isotope_names length ({names}) must match resonance_data length ({isotopes})"
3496 ),
3497 Self::GroupMemberCountMismatch {
3498 group_name,
3499 rd_count,
3500 member_count,
3501 } => write!(
3502 f,
3503 "group '{group_name}': provided {rd_count} ResonanceData but group has {member_count} members"
3504 ),
3505 Self::GroupMemberIsotopeMismatch {
3506 group_name,
3507 expected_z,
3508 expected_a,
3509 got_z,
3510 got_a,
3511 } => write!(
3512 f,
3513 "group '{group_name}': expected Z={expected_z} A={expected_a} but got Z={got_z} A={got_a}"
3514 ),
3515 Self::NonFiniteTemperature(v) => {
3516 write!(f, "temperature must be finite, got {v}")
3517 }
3518 Self::NegativeTemperature(v) => {
3519 write!(f, "temperature must be non-negative, got {v}")
3520 }
3521 Self::InvalidFitEnergyRange(msg) => {
3522 write!(f, "invalid fit_energy_range: {msg}")
3523 }
3524 Self::DensityFreezeBeforeGroups => write!(
3525 f,
3526 "density freeze (with_fix_densities / with_density_free) must be configured \
3527 after with_groups: grouping redefines the density parameters"
3528 ),
3529 }
3530 }
3531}
3532
3533impl std::error::Error for FitConfigError {}
3534
3535#[derive(Debug, Clone)]
3537pub struct SpectrumFitResult {
3538 pub densities: Vec<f64>,
3540 pub uncertainties: Option<Vec<f64>>,
3544 pub reduced_chi_squared: f64,
3546 pub converged: bool,
3548 pub iterations: usize,
3550 pub temperature_k: Option<f64>,
3552 pub temperature_k_unc: Option<f64>,
3564 pub anorm: f64,
3568 pub background: [f64; 3],
3581 pub back_d: Option<f64>,
3589 pub back_f: Option<f64>,
3596 pub t0_us: Option<f64>,
3599 pub l_scale: Option<f64>,
3602 pub energy_scale_flight_path_m: Option<f64>,
3609 pub deviance_per_dof: Option<f64>,
3618 pub baseline: Option<[f64; 3]>,
3635 pub baseline_e_ref_ev: Option<f64>,
3642 pub warnings: Vec<String>,
3651}
3652
3653impl SpectrumFitResult {
3654 pub fn corrected_energies(
3677 &self,
3678 nominal_energies: &[f64],
3679 ) -> Option<Result<Vec<f64>, PipelineError>> {
3680 match (self.t0_us, self.l_scale, self.energy_scale_flight_path_m) {
3681 (Some(t0), Some(l_scale), Some(flight_path_m)) => Some(
3682 nereids_fitting::resolution_calib::corrected_energy_grid(
3683 nominal_energies,
3684 t0,
3685 l_scale,
3686 flight_path_m,
3687 )
3688 .map_err(|e| PipelineError::InvalidParameter(e.to_string())),
3689 ),
3690 _ => None,
3691 }
3692 }
3693}
3694
3695#[cfg(test)]
3696mod tests {
3697 use super::*;
3698 use nereids_endf::resonance::test_support::{
3699 hf178_mlbw_two_resonances, synthetic_single_resonance, u238_single_resonance,
3700 u238_three_resonances,
3701 };
3702 use nereids_fitting::lm::FitModel;
3703 use nereids_fitting::transmission_model::{
3704 EnergyScaleJacobianMethod, EnergyScaleTransmissionModel,
3705 };
3706 use nereids_physics::transmission as phys_transmission;
3707
3708 fn invert_dense(a: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
3712 let n = a.len();
3713 let mut m: Vec<Vec<f64>> = a.to_vec();
3714 let mut inv: Vec<Vec<f64>> = (0..n)
3715 .map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect())
3716 .collect();
3717 for col in 0..n {
3718 let mut piv = col;
3720 for r in (col + 1)..n {
3721 if m[r][col].abs() > m[piv][col].abs() {
3722 piv = r;
3723 }
3724 }
3725 if m[piv][col].abs() < 1e-300 {
3726 return None;
3727 }
3728 m.swap(col, piv);
3729 inv.swap(col, piv);
3730 let d = m[col][col];
3731 for j in 0..n {
3732 m[col][j] /= d;
3733 inv[col][j] /= d;
3734 }
3735 for r in 0..n {
3736 if r == col {
3737 continue;
3738 }
3739 let f = m[r][col];
3740 for j in 0..n {
3741 m[r][j] -= f * m[col][j];
3742 inv[r][j] -= f * inv[col][j];
3743 }
3744 }
3745 }
3746 Some(inv)
3747 }
3748
3749 #[test]
3752 fn detect_transmission_dips_finds_clear_dips() {
3753 let energies: Vec<f64> = (0..120).map(|i| 1.0 + (i as f64) * 0.5).collect();
3754 let mut t = vec![1.0_f64; 120];
3755 for v in t.iter_mut().take(33).skip(28) {
3757 *v = 0.4;
3758 }
3759 for v in t.iter_mut().take(83).skip(78) {
3760 *v = 0.6;
3761 }
3762 let dips = detect_transmission_dips(&t, &energies);
3763 assert_eq!(dips.len(), 2, "expected 2 dips, got {dips:?}");
3764 let mut de: Vec<f64> = dips.iter().map(|&(e, _)| e).collect();
3765 de.sort_by(f64::total_cmp);
3766 assert!(
3768 (de[0] - energies[30]).abs() <= 1.0,
3769 "first dip at {} eV",
3770 de[0]
3771 );
3772 assert!(
3773 (de[1] - energies[80]).abs() <= 1.0,
3774 "second dip at {} eV",
3775 de[1]
3776 );
3777 }
3778
3779 #[test]
3785 fn baseline_reference_energy_honours_fit_energy_range() {
3786 let data = hf178_mlbw_two_resonances();
3787 let mut energies: Vec<f64> = (0..400).map(|i| 8.0 + (i as f64) * 0.0925).collect();
3789 energies.push(3211.0 * 3211.0 / 8.0); let config = UnifiedFitConfig::new(
3791 energies.clone(),
3792 vec![data],
3793 vec!["Hf-178".into()],
3794 293.6,
3795 None,
3796 vec![0.1],
3797 )
3798 .unwrap();
3799 let full = nereids_fitting::transmission_model::baseline_reference_energy(&energies);
3800 assert!((config.baseline_reference_energy() - full).abs() < 1e-6);
3802 let windowed = config.with_fit_energy_range(Some((8.0, 45.0))).unwrap();
3804 let e_ref = windowed.baseline_reference_energy();
3805 assert!(
3806 e_ref > 8.0 && e_ref < 45.0,
3807 "windowed E_ref = {e_ref} eV must lie inside the 8–45 eV fit window"
3808 );
3809 assert!(
3810 (e_ref - full).abs() > 100.0,
3811 "windowed E_ref must differ from the buggy full-grid value {full}"
3812 );
3813 }
3814
3815 #[test]
3820 fn peak_match_energy_scale_seed_identity() {
3821 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
3823 let (t_obs, _sigma) = synthetic_transmission(&data, 0.1, &energies);
3824 let config = UnifiedFitConfig::new(
3825 energies.clone(),
3826 vec![data],
3827 vec!["Hf-178".into()],
3828 293.6,
3829 None,
3830 vec![0.1],
3831 )
3832 .unwrap()
3833 .with_energy_scale(0.0, 1.0, 25.0);
3834 let (t0, l_scale) = peak_match_energy_scale_seed(
3835 &t_obs,
3836 config.energies(),
3837 &config,
3838 25.0,
3839 (-10.0, 10.0),
3840 (0.99, 1.01),
3841 )
3842 .expect("seed should be Some with 2 resonances and detectable dips");
3843 assert!(
3845 t0.abs() < 0.5,
3846 "t0 should be ≈0 for un-shifted data, got {t0}"
3847 );
3848 assert!(
3849 (l_scale - 1.0).abs() < 5e-3,
3850 "L_scale should be ≈1 for un-shifted data, got {l_scale}"
3851 );
3852 }
3853
3854 #[test]
3859 fn peak_match_energy_scale_seed_recovers_nonidentity() {
3860 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
3864 let density = 0.1_f64;
3865 let flight_path = 25.0_f64;
3866 let (t0_true, l_scale_true) = (2.0_f64, 1.006_f64);
3867 let model = EnergyScaleTransmissionModel::new(
3870 Arc::new(vec![data.clone()]),
3871 Arc::new(vec![0]),
3872 Arc::new(vec![1.0]),
3873 293.6,
3874 energies.clone(),
3875 flight_path,
3876 1, 2, None,
3879 );
3880 let t_obs = model.evaluate(&[density, t0_true, l_scale_true]).unwrap();
3881 let config = UnifiedFitConfig::new(
3882 energies.clone(),
3883 vec![data],
3884 vec!["Hf-178".into()],
3885 293.6,
3886 None,
3887 vec![density],
3888 )
3889 .unwrap()
3890 .with_energy_scale(0.0, 1.0, flight_path);
3891 let (t0, l_scale) = peak_match_energy_scale_seed(
3892 &t_obs,
3893 config.energies(),
3894 &config,
3895 flight_path,
3896 (-10.0, 10.0),
3897 (0.99, 1.01),
3898 )
3899 .expect("seed should be Some for a clean shifted two-resonance spectrum");
3900 assert!(
3904 (t0 - t0_true).abs() < 0.6,
3905 "seed should recover t0 ≈ {t0_true}, got {t0}"
3906 );
3907 assert!(
3908 (l_scale - l_scale_true).abs() < 4e-3,
3909 "seed should recover L_scale ≈ {l_scale_true}, got {l_scale}"
3910 );
3911 assert!(
3913 (t0 - t0_true).abs() < (0.0 - t0_true).abs()
3914 && (l_scale - l_scale_true).abs() < (1.0 - l_scale_true).abs(),
3915 "seed must improve on the cold start"
3916 );
3917 }
3918
3919 #[test]
3924 fn peak_match_energy_scale_seed_none_on_featureless_spectrum() {
3925 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
3927 let t_obs: Vec<f64> = energies.iter().map(|&e| (-5.0 / e.sqrt()).exp()).collect();
3929 let config = UnifiedFitConfig::new(
3930 energies.clone(),
3931 vec![data],
3932 vec!["Hf-178".into()],
3933 293.6,
3934 None,
3935 vec![0.1],
3936 )
3937 .unwrap()
3938 .with_energy_scale(0.0, 1.0, 25.0);
3939 assert!(
3940 peak_match_energy_scale_seed(
3941 &t_obs,
3942 config.energies(),
3943 &config,
3944 25.0,
3945 (-10.0, 10.0),
3946 (0.99, 1.01),
3947 )
3948 .is_none(),
3949 "a featureless heavily-absorbing spectrum has no resonance dips ⇒ \
3950 seed must return None (cold-start fallback)"
3951 );
3952 }
3953
3954 #[test]
3961 fn peak_match_energy_scale_seed_handles_duplicate_resonance_energies() {
3962 use nereids_physics::transmission::{SampleParams, forward_model};
3963
3964 let iso_a = synthetic_single_resonance(72, 178, 176.0, 7.8);
3966 let iso_b = synthetic_single_resonance(74, 184, 182.0, 7.8); let iso_c = synthetic_single_resonance(40, 90, 89.0, 16.9); let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
3969 let density = 0.05_f64;
3970 let sample = SampleParams::new(
3971 293.6,
3972 vec![
3973 (iso_a.clone(), density),
3974 (iso_b.clone(), density),
3975 (iso_c.clone(), density),
3976 ],
3977 )
3978 .unwrap();
3979 let t_obs = forward_model(&energies, &sample, None).unwrap();
3980 let config = UnifiedFitConfig::new(
3981 energies.clone(),
3982 vec![iso_a, iso_b, iso_c],
3983 vec!["A".into(), "B".into(), "C".into()],
3984 293.6,
3985 None,
3986 vec![density, density, density],
3987 )
3988 .unwrap()
3989 .with_energy_scale(0.0, 1.0, 25.0);
3990 let (t0, l_scale) = peak_match_energy_scale_seed(
3991 &t_obs,
3992 config.energies(),
3993 &config,
3994 25.0,
3995 (-10.0, 10.0),
3996 (0.99, 1.01),
3997 )
3998 .expect(
3999 "duplicate resonance energies must not collapse match_tol to 0 — \
4000 seed should be Some (was silently None pre-fix)",
4001 );
4002 assert!(t0.abs() < 0.5, "identity calibration: t0 ≈ 0, got {t0}");
4003 assert!(
4004 (l_scale - 1.0).abs() < 5e-3,
4005 "identity calibration: L_scale ≈ 1, got {l_scale}"
4006 );
4007 }
4008
4009 #[test]
4012 fn test_input_data_transmission_n_energies() {
4013 let data = InputData::Transmission {
4014 transmission: vec![0.9, 0.8, 0.7],
4015 uncertainty: vec![0.01, 0.01, 0.01],
4016 };
4017 assert_eq!(data.n_energies(), 3);
4018 assert!(!data.is_counts());
4019 }
4020
4021 #[test]
4022 fn test_input_data_counts_n_energies() {
4023 let data = InputData::Counts {
4024 sample_counts: vec![10.0, 20.0, 30.0, 40.0],
4025 open_beam_counts: vec![100.0, 100.0, 100.0, 100.0],
4026 };
4027 assert_eq!(data.n_energies(), 4);
4028 assert!(data.is_counts());
4029 }
4030
4031 #[test]
4032 fn test_input_data_counts_with_nuisance() {
4033 let data = InputData::CountsWithNuisance {
4034 sample_counts: vec![5.0, 6.0],
4035 flux: vec![100.0, 100.0],
4036 background: vec![0.5, 0.5],
4037 };
4038 assert_eq!(data.n_energies(), 2);
4039 assert!(data.is_counts());
4040 }
4041
4042 #[test]
4043 fn test_solver_config_default_is_auto() {
4044 let cfg = SolverConfig::default();
4045 assert!(matches!(cfg, SolverConfig::Auto));
4046 }
4047
4048 #[test]
4049 fn test_counts_background_config_default() {
4050 let cfg = CountsBackgroundConfig::default();
4051 assert!((cfg.alpha_1_init - 1.0).abs() < f64::EPSILON);
4052 assert!((cfg.alpha_2_init - 1.0).abs() < f64::EPSILON);
4053 assert!(!cfg.fit_alpha_1);
4054 assert!(!cfg.fit_alpha_2);
4055 }
4056
4057 fn synthetic_transmission(
4061 data: &ResonanceData,
4062 true_density: f64,
4063 energies: &[f64],
4064 ) -> (Vec<f64>, Vec<f64>) {
4065 let model = PrecomputedTransmissionModel {
4066 cross_sections: Arc::new(vec![
4067 phys_transmission::broadened_cross_sections(
4068 energies,
4069 std::slice::from_ref(data),
4070 0.0,
4071 None,
4072 None,
4073 )
4074 .unwrap()
4075 .into_iter()
4076 .next()
4077 .unwrap(),
4078 ]),
4079 density_indices: Arc::new(vec![0]),
4080 instrument: None,
4081 resolution_plan: None,
4082 sparse_cubature_plan: None,
4083 sparse_scalar_plan: None,
4084 layout: Arc::new(WorkingGridLayout::identity(energies)),
4085 };
4086 let t = model.evaluate(&[true_density]).unwrap();
4087 let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4088 (t, sigma)
4089 }
4090
4091 fn synthetic_counts(
4093 data: &ResonanceData,
4094 true_density: f64,
4095 energies: &[f64],
4096 i0: f64,
4097 ) -> (Vec<f64>, Vec<f64>) {
4098 let (t, _) = synthetic_transmission(data, true_density, energies);
4099 let open_beam: Vec<f64> = vec![i0; energies.len()];
4100 let sample: Vec<f64> = t.iter().map(|&v| (v * i0).round().max(0.0)).collect();
4101 (sample, open_beam)
4102 }
4103
4104 #[test]
4105 fn test_typed_transmission_lm_recovers_density() {
4106 let data = u238_single_resonance();
4107 let true_density = 0.002;
4108 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4109 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
4110
4111 let config = UnifiedFitConfig::new(
4112 energies,
4113 vec![data],
4114 vec!["U-238".into()],
4115 0.0,
4116 None,
4117 vec![0.001],
4118 )
4119 .unwrap()
4120 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4121
4122 let input = InputData::Transmission {
4123 transmission: t,
4124 uncertainty: sigma,
4125 };
4126
4127 let result = fit_spectrum_typed(&input, &config).unwrap();
4128 assert!(result.converged, "LM should converge");
4129 let fitted = result.densities[0];
4130 assert!(
4131 (fitted - true_density).abs() / true_density < 0.05,
4132 "density: fitted={fitted}, true={true_density}"
4133 );
4134 }
4135
4136 fn table_on_data_grid(config: &UnifiedFitConfig, sigma: Vec<Vec<f64>>) -> PrecomputedXs {
4137 PrecomputedXs {
4138 sigma: Arc::new(sigma),
4139 layout: Arc::new(WorkingGridLayout::identity(config.energies())),
4140 }
4141 }
4142
4143 fn precomputed_xs_fixture() -> (UnifiedFitConfig, InputData) {
4146 let data = u238_single_resonance();
4147 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
4148 let (t, sigma) = synthetic_transmission(&data, 0.001, &energies);
4149 let config = UnifiedFitConfig::new(
4150 energies,
4151 vec![data],
4152 vec!["U-238".into()],
4153 0.0,
4154 None,
4155 vec![0.001],
4156 )
4157 .unwrap()
4158 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4159 let input = InputData::Transmission {
4160 transmission: t,
4161 uncertainty: sigma,
4162 };
4163 (config, input)
4164 }
4165
4166 #[test]
4167 fn test_fit_rejects_empty_precomputed_cross_sections() {
4168 let (config, input) = precomputed_xs_fixture();
4171 let config = config
4172 .clone()
4173 .with_precomputed_cross_sections(table_on_data_grid(&config, Vec::new()));
4174 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4175 assert!(
4176 matches!(err, PipelineError::ShapeMismatch(_)),
4177 "expected ShapeMismatch, got {err:?}"
4178 );
4179 assert!(err.to_string().contains("must not be empty"));
4180 }
4181
4182 #[test]
4183 fn test_fit_rejects_wrong_row_count_precomputed_cross_sections() {
4184 let (config, input) = precomputed_xs_fixture();
4186 let n_e = config.energies().len();
4187 let bad = vec![vec![1.0; n_e], vec![1.0; n_e]];
4188 let config = config
4189 .clone()
4190 .with_precomputed_cross_sections(table_on_data_grid(&config, bad));
4191 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4192 assert!(
4193 matches!(err, PipelineError::ShapeMismatch(_)),
4194 "expected ShapeMismatch, got {err:?}"
4195 );
4196 assert!(err.to_string().contains("rows"));
4197 }
4198
4199 #[test]
4200 fn test_fit_rejects_wrong_energy_length_precomputed_cross_sections() {
4201 let (config, input) = precomputed_xs_fixture();
4203 let n_e = config.energies().len();
4204 let bad = vec![vec![1.0; n_e + 3]];
4205 let config = config
4206 .clone()
4207 .with_precomputed_cross_sections(table_on_data_grid(&config, bad));
4208 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4209 assert!(
4210 matches!(err, PipelineError::ShapeMismatch(_)),
4211 "expected ShapeMismatch, got {err:?}"
4212 );
4213 assert!(err.to_string().contains("its grid has"));
4214 }
4215
4216 #[test]
4217 fn test_fit_accepts_correct_precomputed_cross_sections() {
4218 let (config, input) = precomputed_xs_fixture();
4221 let n_e = config.energies().len();
4222 let xs = phys_transmission::broadened_cross_sections(
4223 config.energies(),
4224 config.resonance_data(),
4225 0.0,
4226 None,
4227 None,
4228 )
4229 .unwrap();
4230 assert_eq!(xs.len(), 1);
4231 assert_eq!(xs[0].len(), n_e);
4232 let config = config
4233 .clone()
4234 .with_precomputed_cross_sections(table_on_data_grid(&config, xs));
4235 let result = fit_spectrum_typed(&input, &config);
4238 assert!(
4239 result.is_ok(),
4240 "correctly-shaped precomputed XS must pass validation, got {result:?}"
4241 );
4242 }
4243
4244 #[test]
4245 fn test_fit_rejects_non_finite_precomputed_cross_sections() {
4246 let (config, input) = precomputed_xs_fixture();
4252 let n_e = config.energies().len();
4253 let bad = vec![vec![f64::NAN; n_e]];
4254 let config = config
4255 .clone()
4256 .with_precomputed_cross_sections(table_on_data_grid(&config, bad));
4257 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4258 assert!(
4259 matches!(err, PipelineError::ShapeMismatch(_)),
4260 "expected ShapeMismatch, got {err:?}"
4261 );
4262 assert!(
4263 err.to_string().contains("non-finite"),
4264 "error should mention non-finite σ, got: {err}"
4265 );
4266
4267 let (config, input) = precomputed_xs_fixture();
4269 let mut row = vec![1.0; n_e];
4270 row[n_e / 2] = f64::INFINITY;
4271 let config = config
4272 .clone()
4273 .with_precomputed_cross_sections(table_on_data_grid(&config, vec![row]));
4274 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4275 assert!(
4276 matches!(err, PipelineError::ShapeMismatch(_)),
4277 "expected ShapeMismatch for +inf σ, got {err:?}"
4278 );
4279 }
4280
4281 #[test]
4282 fn test_evaluate_jacobian_rejects_malformed_precomputed_cross_sections() {
4283 let (config, _input) = precomputed_xs_fixture();
4288 let n_e = config.energies().len();
4289 let flux = vec![1.0; n_e];
4290 let background = vec![0.0; n_e];
4291
4292 let empty = config
4295 .clone()
4296 .with_precomputed_cross_sections(table_on_data_grid(&config, Vec::new()));
4297 match evaluate_jacobian_and_fisher(&empty, &flux, &background) {
4298 Err(PipelineError::ShapeMismatch(msg)) => {
4299 assert!(msg.contains("must not be empty"), "got: {msg}");
4300 }
4301 Err(other) => panic!("expected ShapeMismatch for empty XS, got {other:?}"),
4302 Ok(_) => panic!("empty precomputed XS must be rejected"),
4303 }
4304
4305 let nan = config
4307 .clone()
4308 .with_precomputed_cross_sections(table_on_data_grid(
4309 &config,
4310 vec![vec![f64::NAN; n_e]],
4311 ));
4312 match evaluate_jacobian_and_fisher(&nan, &flux, &background) {
4313 Err(PipelineError::ShapeMismatch(msg)) => {
4314 assert!(msg.contains("non-finite"), "got: {msg}");
4315 }
4316 Err(other) => panic!("expected ShapeMismatch for NaN σ, got {other:?}"),
4317 Ok(_) => panic!("non-finite precomputed σ must be rejected"),
4318 }
4319 }
4320
4321 #[test]
4322 fn test_extract_result_drops_uncertainties_when_unconverged() {
4323 let data = u238_single_resonance();
4324 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
4325 let config = UnifiedFitConfig::new(
4326 energies,
4327 vec![data],
4328 vec!["U-238".into()],
4329 293.6,
4330 None,
4331 vec![0.001],
4332 )
4333 .unwrap();
4334
4335 let result = LmResult {
4336 chi_squared: 1.0,
4337 reduced_chi_squared: 1.0,
4338 iterations: 5,
4339 converged: false,
4340 params: vec![0.001],
4341 covariance: Some(lm::FlatMatrix::zeros(1, 1)),
4342 uncertainties: Some(vec![0.123]),
4343 };
4344
4345 let extracted = extract_result(&config, &result, 1, &[0], None, None).unwrap();
4347 assert!(!extracted.converged);
4348 assert!(
4349 extracted.uncertainties.is_none(),
4350 "pipeline must not surface uncertainties from an unconverged fit"
4351 );
4352 assert!(extracted.temperature_k_unc.is_none());
4353 }
4354
4355 #[test]
4356 fn test_typed_counts_kl_recovers_density() {
4357 let data = u238_single_resonance();
4358 let true_density = 0.002;
4359 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4360 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4361
4362 let config = UnifiedFitConfig::new(
4363 energies,
4364 vec![data],
4365 vec!["U-238".into()],
4366 0.0,
4367 None,
4368 vec![0.001],
4369 )
4370 .unwrap()
4371 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4372
4373 let input = InputData::Counts {
4374 sample_counts: sample,
4375 open_beam_counts: open_beam,
4376 };
4377
4378 let result = fit_spectrum_typed(&input, &config).unwrap();
4379 assert!(result.converged, "KL on counts should converge");
4380 let fitted = result.densities[0];
4381 assert!(
4382 (fitted - true_density).abs() / true_density < 0.10,
4383 "density: fitted={fitted}, true={true_density}"
4384 );
4385 }
4386
4387 #[test]
4388 fn test_typed_counts_kl_low_counts_recovers_density() {
4389 let data = u238_single_resonance();
4391 let true_density = 0.0005;
4392 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4393 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 10.0);
4394
4395 let config = UnifiedFitConfig::new(
4396 energies,
4397 vec![data],
4398 vec!["U-238".into()],
4399 0.0,
4400 None,
4401 vec![0.001],
4402 )
4403 .unwrap()
4404 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4405
4406 let input = InputData::Counts {
4407 sample_counts: sample,
4408 open_beam_counts: open_beam,
4409 };
4410
4411 let result = fit_spectrum_typed(&input, &config).unwrap();
4412 assert!(result.converged, "KL on low counts should converge");
4413 let fitted = result.densities[0];
4414 assert!(
4416 (fitted - true_density).abs() / true_density < 0.30,
4417 "density: fitted={fitted}, true={true_density}"
4418 );
4419 }
4420
4421 #[test]
4422 fn test_typed_transmission_kl_is_rejected() {
4423 let data = u238_single_resonance();
4424 let true_density = 0.0005;
4425 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4426 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
4427
4428 let config = UnifiedFitConfig::new(
4429 energies,
4430 vec![data],
4431 vec!["U-238".into()],
4432 0.0,
4433 None,
4434 vec![0.001],
4435 )
4436 .unwrap()
4437 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4438
4439 let input = InputData::Transmission {
4440 transmission: t,
4441 uncertainty: sigma,
4442 };
4443
4444 let error = fit_spectrum_typed(&input, &config)
4445 .expect_err("fractional transmission must not use a Poisson count objective");
4446 let message = error.to_string();
4447 assert!(
4448 message.contains("transmission")
4449 && message.contains("Poisson")
4450 && message.contains("least-squares"),
4451 "rejection must explain the valid transmission engine, got: {message}"
4452 );
4453 }
4454
4455 #[test]
4456 fn test_typed_counts_lm_is_rejected() {
4457 let data = u238_single_resonance();
4460 let true_density = 0.0005;
4461 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4462 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4463
4464 let config = UnifiedFitConfig::new(
4465 energies,
4466 vec![data],
4467 vec!["U-238".into()],
4468 0.0,
4469 None,
4470 vec![0.001],
4471 )
4472 .unwrap()
4473 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4474
4475 let input = InputData::Counts {
4476 sample_counts: sample,
4477 open_beam_counts: open_beam,
4478 };
4479
4480 let error = fit_spectrum_typed(&input, &config)
4481 .expect_err("raw counts must not be silently converted to transmission");
4482 let message = error.to_string();
4483 assert!(
4484 message.contains("counts")
4485 && message.contains("least-squares")
4486 && message.contains("Poisson"),
4487 "rejection must explain the valid counts engine, got: {message}"
4488 );
4489 }
4490
4491 #[test]
4492 fn transmission_rejects_counts_background_config() {
4493 let data = u238_single_resonance();
4494 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
4495 let (transmission, uncertainty) = synthetic_transmission(&data, 0.0005, &energies);
4496 let config = UnifiedFitConfig::new(
4497 energies,
4498 vec![data],
4499 vec!["U-238".into()],
4500 0.0,
4501 None,
4502 vec![0.001],
4503 )
4504 .unwrap()
4505 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
4506 .with_counts_background(CountsBackgroundConfig::default());
4507 let input = InputData::Transmission {
4508 transmission,
4509 uncertainty,
4510 };
4511
4512 let error = fit_spectrum_typed(&input, &config)
4513 .expect_err("a counts-background request must not be ignored on transmission");
4514 let message = error.to_string();
4515 assert!(
4516 message.contains("counts background") && message.contains("transmission"),
4517 "rejection must name the domain mismatch, got: {message}"
4518 );
4519 }
4520
4521 #[test]
4522 fn test_typed_auto_solver_selects_kl_for_counts() {
4523 let data = u238_single_resonance();
4524 let true_density = 0.0005;
4525 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4526 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4527
4528 let config = UnifiedFitConfig::new(
4530 energies,
4531 vec![data],
4532 vec!["U-238".into()],
4533 0.0,
4534 None,
4535 vec![0.001],
4536 )
4537 .unwrap(); let input = InputData::Counts {
4540 sample_counts: sample,
4541 open_beam_counts: open_beam,
4542 };
4543
4544 let result = fit_spectrum_typed(&input, &config).unwrap();
4545 assert!(
4546 result.converged,
4547 "Auto solver on counts should use KL and converge"
4548 );
4549 }
4550
4551 #[test]
4552 fn test_typed_transmission_with_background_lm() {
4553 let data = u238_single_resonance();
4554 let true_density = 0.0005;
4555 let true_anorm = 0.95;
4556 let true_back_a = 0.02;
4557 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4558
4559 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
4561 let t_bg: Vec<f64> = t_pure
4562 .iter()
4563 .map(|&v| true_anorm * v + true_back_a)
4564 .collect();
4565 let sigma: Vec<f64> = t_bg.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4566
4567 let config = UnifiedFitConfig::new(
4568 energies,
4569 vec![data],
4570 vec!["U-238".into()],
4571 0.0,
4572 None,
4573 vec![0.001],
4574 )
4575 .unwrap()
4576 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
4577 max_iter: 500,
4578 ..LmConfig::default()
4579 }))
4580 .with_transmission_background(BackgroundConfig::default());
4581
4582 let input = InputData::Transmission {
4583 transmission: t_bg,
4584 uncertainty: sigma,
4585 };
4586
4587 let result = fit_spectrum_typed(&input, &config).unwrap();
4588 assert!(
4589 result.converged,
4590 "LM+BG should converge on noiseless synthetic data (chi2r={}, iter={})",
4591 result.reduced_chi_squared, result.iterations
4592 );
4593 assert!(
4594 (result.densities[0] - true_density).abs() / true_density < 0.05,
4595 "density: fitted={}, true={true_density}",
4596 result.densities[0]
4597 );
4598 assert!(
4599 (result.anorm - true_anorm).abs() / true_anorm < 0.05,
4600 "anorm: fitted={}, true={true_anorm}",
4601 result.anorm
4602 );
4603 }
4604
4605 #[test]
4606 fn test_typed_counts_with_nuisance_rejects_lm() {
4607 let data = u238_single_resonance();
4608 let energies: Vec<f64> = (0..10).map(|i| 1.0 + (i as f64) * 0.5).collect();
4609
4610 let config = UnifiedFitConfig::new(
4611 energies,
4612 vec![data],
4613 vec!["U-238".into()],
4614 0.0,
4615 None,
4616 vec![0.001],
4617 )
4618 .unwrap()
4619 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4620
4621 let input = InputData::CountsWithNuisance {
4622 sample_counts: vec![10.0; 10],
4623 flux: vec![100.0; 10],
4624 background: vec![0.0; 10],
4625 };
4626
4627 let result = fit_spectrum_typed(&input, &config);
4628 assert!(result.is_err(), "CountsWithNuisance + LM should error");
4629 }
4630
4631 fn seeded_gaussian(n: usize, seed: u64) -> Vec<f64> {
4635 let mut state = seed | 1;
4636 let mut unif = || {
4637 state = state
4638 .wrapping_mul(6364136223846793005)
4639 .wrapping_add(1442695040888963407);
4640 (((state >> 11) as f64) / ((1u64 << 53) as f64)).clamp(1e-12, 1.0 - 1e-12)
4641 };
4642 (0..n)
4643 .map(|_| {
4644 let u1 = unif();
4645 let u2 = unif();
4646 (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos()
4647 })
4648 .collect()
4649 }
4650
4651 fn synthetic_transmission_at_temp(
4652 data: &ResonanceData,
4653 true_density: f64,
4654 temperature_k: f64,
4655 energies: &[f64],
4656 ) -> (Vec<f64>, Vec<f64>) {
4657 let xs = phys_transmission::broadened_cross_sections(
4658 energies,
4659 std::slice::from_ref(data),
4660 temperature_k,
4661 None,
4662 None,
4663 )
4664 .unwrap();
4665 let model = PrecomputedTransmissionModel {
4666 cross_sections: Arc::new(xs),
4667 density_indices: Arc::new(vec![0]),
4668 instrument: None,
4669 resolution_plan: None,
4670 sparse_cubature_plan: None,
4671 sparse_scalar_plan: None,
4672 layout: Arc::new(WorkingGridLayout::identity(energies)),
4673 };
4674 let t = model.evaluate(&[true_density]).unwrap();
4675 let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4676 (t, sigma)
4677 }
4678
4679 #[test]
4680 fn test_typed_lm_with_temperature() {
4681 let data = u238_single_resonance();
4682 let true_density = 0.0005;
4683 let true_temp = 350.0;
4684 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4685 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
4686
4687 let config = UnifiedFitConfig::new(
4688 energies,
4689 vec![data],
4690 vec!["U-238".into()],
4691 300.0, None,
4693 vec![0.001],
4694 )
4695 .unwrap()
4696 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
4697 max_iter: 500,
4698 ..LmConfig::default()
4699 }))
4700 .with_fit_temperature(true);
4701
4702 let input = InputData::Transmission {
4703 transmission: t,
4704 uncertainty: sigma,
4705 };
4706
4707 let result = fit_spectrum_typed(&input, &config).unwrap();
4708
4709 let fitted_density = result.densities[0];
4711 assert!(
4712 (fitted_density - true_density).abs() / true_density < 0.01,
4713 "density: fitted={fitted_density}, true={true_density}, ratio={}",
4714 (fitted_density - true_density).abs() / true_density,
4715 );
4716
4717 let fitted_temp = result
4719 .temperature_k
4720 .expect("temperature_k should be Some when fit_temperature=true");
4721 assert!(
4722 (fitted_temp - true_temp).abs() < 1.0,
4723 "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
4724 (fitted_temp - true_temp).abs(),
4725 );
4726 }
4727
4728 #[test]
4729 fn test_typed_lm_with_temperature_and_background() {
4730 let data = u238_single_resonance();
4731 let true_density = 0.0005;
4732 let true_temp = 350.0;
4733 let true_b0 = 0.012;
4734 let true_b1 = 0.008;
4735 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4736 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
4737 let measured_t: Vec<f64> = t
4738 .iter()
4739 .zip(energies.iter())
4740 .map(|(&ti, &e)| ti + true_b0 + true_b1 / e.sqrt())
4741 .collect();
4742
4743 let config = UnifiedFitConfig::new(
4744 energies,
4745 vec![data],
4746 vec!["U-238".into()],
4747 300.0,
4748 None,
4749 vec![0.001],
4750 )
4751 .unwrap()
4752 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
4753 max_iter: 120,
4754 ..LmConfig::default()
4755 }))
4756 .with_fit_temperature(true)
4757 .with_transmission_background(BackgroundConfig::default());
4758
4759 let input = InputData::Transmission {
4760 transmission: measured_t,
4761 uncertainty: sigma,
4762 };
4763
4764 let result = fit_spectrum_typed(&input, &config).unwrap();
4765
4766 assert!(result.converged, "fit did not converge: {result:?}");
4767 assert!(
4768 result.iterations <= 80,
4769 "expected LM background+temperature fit to converge well before max_iter; got {}",
4770 result.iterations,
4771 );
4772
4773 let fitted_density = result.densities[0];
4774 assert!(
4775 (fitted_density - true_density).abs() / true_density < 0.02,
4776 "density: fitted={fitted_density}, true={true_density}, ratio={}",
4777 (fitted_density - true_density).abs() / true_density,
4778 );
4779
4780 let fitted_temp = result
4781 .temperature_k
4782 .expect("temperature_k should be Some when fit_temperature=true");
4783 assert!(
4784 (fitted_temp - true_temp).abs() < 3.0,
4785 "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
4786 (fitted_temp - true_temp).abs(),
4787 );
4788
4789 assert!(
4790 (result.background[0] - true_b0).abs() < 5e-3,
4791 "background b0: fitted={}, true={}",
4792 result.background[0],
4793 true_b0,
4794 );
4795 assert!(
4796 (result.background[1] - true_b1).abs() < 5e-3,
4797 "background b1: fitted={}, true={}",
4798 result.background[1],
4799 true_b1,
4800 );
4801 }
4802
4803 #[test]
4807 fn test_grouped_fit_spectrum_round_trip() {
4808 use nereids_core::types::IsotopeGroup;
4809
4810 let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
4812 let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);
4813
4814 let iso1 = nereids_core::types::Isotope::new(92, 235).unwrap();
4816 let iso2 = nereids_core::types::Isotope::new(92, 238).unwrap();
4817 let group =
4818 IsotopeGroup::custom("U (60/40)".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();
4819
4820 let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
4821 let true_density = 0.0005;
4822
4823 let sample = nereids_physics::transmission::SampleParams::new(
4825 0.0,
4826 vec![
4827 (rd1.clone(), true_density * 0.6),
4828 (rd2.clone(), true_density * 0.4),
4829 ],
4830 )
4831 .unwrap();
4832 let transmission =
4833 nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
4834 let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();
4835
4836 let config = UnifiedFitConfig::new(
4838 energies.clone(),
4839 vec![rd1.clone()],
4840 vec!["placeholder".into()],
4841 0.0,
4842 None,
4843 vec![0.001],
4844 )
4845 .unwrap()
4846 .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
4847 .unwrap()
4848 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4849
4850 let input = InputData::Transmission {
4851 transmission,
4852 uncertainty,
4853 };
4854
4855 let result = fit_spectrum_typed(&input, &config).unwrap();
4856
4857 assert_eq!(result.densities.len(), 1, "should have 1 group density");
4859 let fitted = result.densities[0];
4860 let rel_error = (fitted - true_density).abs() / true_density;
4861 assert!(
4862 rel_error < 0.01,
4863 "group density: fitted={fitted}, true={true_density}, rel_error={rel_error}"
4864 );
4865 assert!(result.converged, "fit should converge");
4866 }
4867
4868 #[test]
4869 fn test_grouped_lm_with_temperature_and_background_noiseless() {
4870 use nereids_core::types::IsotopeGroup;
4871
4872 let rd1 = synthetic_single_resonance(72, 176, 8.5, 5.0);
4873 let rd2 = synthetic_single_resonance(72, 178, 17.0, 7.5);
4874 let rd3 = synthetic_single_resonance(72, 180, 29.0, 6.0);
4875
4876 let hf176 = nereids_core::types::Isotope::new(72, 176).unwrap();
4877 let hf178 = nereids_core::types::Isotope::new(72, 178).unwrap();
4878 let hf180 = nereids_core::types::Isotope::new(72, 180).unwrap();
4879 let group = IsotopeGroup::custom(
4880 "Hf-like (3 member)".into(),
4881 vec![(hf176, 0.2), (hf178, 0.5), (hf180, 0.3)],
4882 )
4883 .unwrap();
4884
4885 let energies: Vec<f64> = (0..300).map(|i| 1.0 + (49.0 * i as f64) / 299.0).collect();
4886 let true_density = 0.001;
4887 let true_temp = 400.0;
4888 let true_b0 = 0.012;
4889 let true_b1 = 0.008;
4890
4891 let sample = nereids_physics::transmission::SampleParams::new(
4892 true_temp,
4893 vec![
4894 (rd1.clone(), true_density * 0.2),
4895 (rd2.clone(), true_density * 0.5),
4896 (rd3.clone(), true_density * 0.3),
4897 ],
4898 )
4899 .unwrap();
4900 let pure_t =
4901 nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
4902 let measured_t: Vec<f64> = pure_t
4903 .iter()
4904 .zip(energies.iter())
4905 .map(|(&t, &e)| t + true_b0 + true_b1 / e.sqrt())
4906 .collect();
4907 let sigma = vec![0.001; energies.len()];
4908
4909 let config = UnifiedFitConfig::new(
4910 energies.clone(),
4911 vec![rd1.clone()],
4912 vec!["placeholder".into()],
4913 293.6,
4914 None,
4915 vec![0.0008],
4916 )
4917 .unwrap()
4918 .with_groups(&[(&group, &[rd1, rd2, rd3])], vec![0.0008])
4919 .unwrap()
4920 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
4921 max_iter: 200,
4922 ..LmConfig::default()
4923 }))
4924 .with_fit_temperature(true)
4925 .with_transmission_background(BackgroundConfig::default());
4926
4927 let input = InputData::Transmission {
4928 transmission: measured_t,
4929 uncertainty: sigma,
4930 };
4931
4932 let result = fit_spectrum_typed(&input, &config).unwrap();
4933
4934 assert!(result.converged, "fit did not converge: {result:?}");
4935 assert!(
4938 (result.densities[0] - true_density).abs() / true_density < 0.01,
4939 "density: fitted={}, true={true_density}",
4940 result.densities[0]
4941 );
4942 let fitted_temp = result
4943 .temperature_k
4944 .expect("temperature_k should be Some when fit_temperature=true");
4945 assert!(
4946 (fitted_temp - true_temp).abs() < 8.0,
4947 "temperature: fitted={fitted_temp}, true={true_temp}",
4948 );
4949 let e_mid: f64 = 10.0;
4953 let bg_total = (result.anorm - 1.0)
4954 + result.background[0]
4955 + result.background[1] / e_mid.sqrt()
4956 + result.background[2] * e_mid.sqrt();
4957 let true_bg_mid = true_b0 + true_b1 / e_mid.sqrt();
4958 assert!(
4959 (bg_total - true_bg_mid).abs() < 0.02,
4960 "total bg at E={e_mid}: fitted={bg_total:.6}, true={true_bg_mid:.6}",
4961 );
4962 }
4963
4964 #[test]
4967 fn test_kl_counts_returns_density_uncertainty() {
4968 let data = u238_single_resonance();
4969 let true_density = 0.002;
4970 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4971 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4972
4973 let config = UnifiedFitConfig::new(
4974 energies,
4975 vec![data],
4976 vec!["U-238".into()],
4977 0.0,
4978 None,
4979 vec![0.001],
4980 )
4981 .unwrap()
4982 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4983
4984 let input = InputData::Counts {
4985 sample_counts: sample,
4986 open_beam_counts: open_beam,
4987 };
4988 let result = fit_spectrum_typed(&input, &config).unwrap();
4989 assert!(result.converged);
4990 let unc = result
4991 .uncertainties
4992 .as_ref()
4993 .expect("KL 1D fit should return density uncertainties");
4994 assert_eq!(unc.len(), 1);
4995 assert!(
4996 unc[0].is_finite() && unc[0] > 0.0,
4997 "density unc = {}",
4998 unc[0]
4999 );
5000 assert!(
5001 unc[0] < result.densities[0],
5002 "unc ({}) should be < density ({}) for high-count data",
5003 unc[0],
5004 result.densities[0]
5005 );
5006 }
5007
5008 #[test]
5009 fn test_kl_counts_returns_temperature_uncertainty() {
5010 let data = u238_single_resonance();
5011 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5012 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5013
5014 let config = UnifiedFitConfig::new(
5015 energies,
5016 vec![data],
5017 vec!["U-238".into()],
5018 350.0,
5019 None,
5020 vec![0.0005],
5021 )
5022 .unwrap()
5023 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5024 .with_fit_temperature(true);
5025
5026 let input = InputData::Counts {
5027 sample_counts: sample,
5028 open_beam_counts: open_beam,
5029 };
5030 let result = fit_spectrum_typed(&input, &config).unwrap();
5031 assert!(result.converged);
5032 let unc = result
5033 .uncertainties
5034 .as_ref()
5035 .expect("KL+temp fit should return density uncertainties");
5036 assert!(
5037 unc[0].is_finite() && unc[0] > 0.0,
5038 "density unc = {}",
5039 unc[0]
5040 );
5041 let t_unc = result
5042 .temperature_k_unc
5043 .expect("KL+temp fit should return temperature uncertainty");
5044 assert!(
5045 t_unc.is_finite() && t_unc > 0.0,
5046 "temperature unc = {t_unc}"
5047 );
5048 }
5049
5050 #[test]
5053 fn test_scale_by_chi2_inflates_joint_poisson_uncertainty() {
5054 let data = u238_single_resonance();
5055 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5056 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5057 let input = InputData::Counts {
5058 sample_counts: sample,
5059 open_beam_counts: open_beam,
5060 };
5061
5062 let make_config = |scale: bool| {
5063 UnifiedFitConfig::new(
5064 energies.clone(),
5065 vec![data.clone()],
5066 vec!["U-238".into()],
5067 350.0,
5068 None,
5069 vec![0.0005],
5070 )
5071 .unwrap()
5072 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5073 .with_fit_temperature(true)
5074 .with_scale_by_chi2(scale)
5075 };
5076
5077 let unscaled = fit_spectrum_typed(&input, &make_config(false)).unwrap();
5078 let scaled = fit_spectrum_typed(&input, &make_config(true)).unwrap();
5079 assert!(unscaled.converged && scaled.converged);
5080
5081 let dpd = unscaled
5084 .deviance_per_dof
5085 .expect("joint-Poisson must report deviance_per_dof");
5086 assert!(dpd.is_finite() && dpd > 0.0, "deviance_per_dof = {dpd}");
5087 let factor = dpd.sqrt();
5088 assert!(
5093 (factor - 1.0).abs() > 0.01,
5094 "expected D/dof to move σ by >1%, got factor {factor}"
5095 );
5096
5097 let t_unscaled = unscaled.temperature_k_unc.expect("σ_T unscaled");
5099 let t_scaled = scaled.temperature_k_unc.expect("σ_T scaled");
5100 let rel_t = (t_scaled - t_unscaled * factor).abs() / (t_unscaled * factor);
5101 assert!(
5102 rel_t < 1e-6,
5103 "σ_T: scaled {t_scaled} must equal unscaled {t_unscaled} × {factor}"
5104 );
5105
5106 let d_unscaled = unscaled.uncertainties.as_ref().expect("density σ unscaled")[0];
5108 let d_scaled = scaled.uncertainties.as_ref().expect("density σ scaled")[0];
5109 let rel_d = (d_scaled - d_unscaled * factor).abs() / (d_unscaled * factor);
5110 assert!(
5111 rel_d < 1e-6,
5112 "density σ: scaled {d_scaled} must equal unscaled {d_unscaled} × {factor}"
5113 );
5114
5115 let default_cfg = UnifiedFitConfig::new(
5118 energies.clone(),
5119 vec![data.clone()],
5120 vec!["U-238".into()],
5121 350.0,
5122 None,
5123 vec![0.0005],
5124 )
5125 .unwrap()
5126 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5127 .with_fit_temperature(true);
5128 let default_run = fit_spectrum_typed(&input, &default_cfg).unwrap();
5129 assert_eq!(
5130 default_run.temperature_k_unc, unscaled.temperature_k_unc,
5131 "default (flag absent) must equal scale_by_chi2=false"
5132 );
5133 }
5134
5135 #[test]
5143 fn test_scale_by_chi2_joint_poisson_underfit_grows_sigma() {
5144 let data = u238_single_resonance();
5145 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5146 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5147 let sample: Vec<f64> = sample
5152 .iter()
5153 .enumerate()
5154 .map(|(i, &c)| {
5155 let k = if i % 2 == 0 { 1.30 } else { 0.70 };
5156 (c * k).round().max(0.0)
5157 })
5158 .collect();
5159 let input = InputData::Counts {
5160 sample_counts: sample,
5161 open_beam_counts: open_beam,
5162 };
5163 let make_config = |scale: bool| {
5164 UnifiedFitConfig::new(
5165 energies.clone(),
5166 vec![data.clone()],
5167 vec!["U-238".into()],
5168 350.0,
5169 None,
5170 vec![0.0005],
5171 )
5172 .unwrap()
5173 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5174 .with_fit_temperature(true)
5175 .with_scale_by_chi2(scale)
5176 };
5177 let unscaled = fit_spectrum_typed(&input, &make_config(false)).unwrap();
5178 let scaled = fit_spectrum_typed(&input, &make_config(true)).unwrap();
5179 assert!(unscaled.converged && scaled.converged);
5180 let dpd = unscaled
5181 .deviance_per_dof
5182 .expect("joint-Poisson must report deviance_per_dof");
5183 assert!(
5184 dpd > 1.0,
5185 "zig-zag counts must under-fit (D/dof > 1), got {dpd}"
5186 );
5187 let t_unscaled = unscaled.temperature_k_unc.expect("σ_T unscaled");
5188 let t_scaled = scaled.temperature_k_unc.expect("σ_T scaled");
5189 assert!(
5190 t_scaled > t_unscaled,
5191 "under-fit (D/dof {dpd} > 1) must GROW σ_T: scaled {t_scaled} \
5192 vs unscaled {t_unscaled}"
5193 );
5194 }
5195
5196 #[test]
5197 fn test_kl_counts_with_background_returns_uncertainty() {
5198 let data = u238_single_resonance();
5199 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5200 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5201
5202 let config = UnifiedFitConfig::new(
5203 energies,
5204 vec![data],
5205 vec!["U-238".into()],
5206 300.0,
5207 None,
5208 vec![0.0005],
5209 )
5210 .unwrap()
5211 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5212 .with_fit_temperature(true)
5213 .with_transmission_background(BackgroundConfig::default());
5214
5215 let input = InputData::Counts {
5216 sample_counts: sample,
5217 open_beam_counts: open_beam,
5218 };
5219 let result = fit_spectrum_typed(&input, &config).unwrap();
5220 assert!(result.converged);
5221 let unc = result
5222 .uncertainties
5223 .as_ref()
5224 .expect("KL+bg fit should return density uncertainties");
5225 assert!(
5226 unc[0].is_finite() && unc[0] > 0.0,
5227 "density unc = {}",
5228 unc[0]
5229 );
5230 }
5231
5232 #[test]
5233 fn test_lm_uncertainty_not_regressed() {
5234 let data = u238_single_resonance();
5235 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5236 let (t_clean, sigma) = synthetic_transmission(&data, 0.001, &energies);
5237 let t: Vec<f64> = t_clean
5244 .iter()
5245 .enumerate()
5246 .map(|(i, &v)| v * (1.0 + 0.002 * (7.3 * i as f64).sin()))
5247 .collect();
5248
5249 let config = UnifiedFitConfig::new(
5250 energies,
5251 vec![data],
5252 vec!["U-238".into()],
5253 0.0,
5254 None,
5255 vec![0.0005],
5256 )
5257 .unwrap()
5258 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
5259
5260 let input = InputData::Transmission {
5261 transmission: t,
5262 uncertainty: sigma,
5263 };
5264 let result = fit_spectrum_typed(&input, &config).unwrap();
5265 assert!(result.converged);
5266 let unc = result
5267 .uncertainties
5268 .as_ref()
5269 .expect("LM should still return uncertainties");
5270 assert!(unc[0].is_finite() && unc[0] > 0.0);
5271 }
5272
5273 #[test]
5291 fn test_energy_scale_with_temperature_recovers_all_three() {
5292 let data = u238_three_resonances();
5295 let flight_path = 25.0_f64;
5296 let true_density = 0.002;
5297 let true_temp = 340.0;
5298 let true_t0 = 0.6_f64; let true_l_scale = 1.004_f64;
5300 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5302
5303 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5305 &nominal,
5306 true_t0,
5307 true_l_scale,
5308 flight_path,
5309 )
5310 .unwrap();
5311 let (t_obs, sigma) =
5313 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5314
5315 let config = UnifiedFitConfig::new(
5316 nominal,
5317 vec![data],
5318 vec!["U-238".into()],
5319 true_temp - 40.0, None,
5321 vec![true_density],
5322 )
5323 .unwrap()
5324 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5325 .with_fit_temperature(true)
5326 .with_energy_scale(0.0, 1.0, flight_path); let input = InputData::Transmission {
5329 transmission: t_obs,
5330 uncertainty: sigma,
5331 };
5332 let result = fit_spectrum_typed(&input, &config).expect("joint fit runs");
5333 assert!(result.converged, "joint (t0,L_scale,T) fit should converge");
5334
5335 let t0 = result.t0_us.expect("t0_us populated");
5336 let ls = result.l_scale.expect("l_scale populated");
5337 let temp = result.temperature_k.expect("temperature_k populated");
5338 assert!(
5339 (t0 - true_t0).abs() < 0.05,
5340 "t0: fitted={t0}, true={true_t0}"
5341 );
5342 assert!(
5343 (ls - true_l_scale).abs() / true_l_scale < 1e-3,
5344 "l_scale: fitted={ls}, true={true_l_scale}"
5345 );
5346 assert!(
5347 (temp - true_temp).abs() < 3.0,
5348 "temperature: fitted={temp}, true={true_temp}"
5349 );
5350 }
5351
5352 #[test]
5355 fn test_energy_scale_with_temperature_does_not_enable_transmission_kl() {
5356 let data = u238_three_resonances();
5357 let flight_path = 25.0_f64;
5358 let true_density = 0.002;
5359 let true_temp = 340.0;
5360 let true_t0 = 0.6_f64;
5361 let true_l_scale = 1.004_f64;
5362 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5363 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5364 &nominal,
5365 true_t0,
5366 true_l_scale,
5367 flight_path,
5368 )
5369 .unwrap();
5370 let (t_obs, sigma) =
5371 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5372
5373 let config = UnifiedFitConfig::new(
5374 nominal,
5375 vec![data],
5376 vec!["U-238".into()],
5377 true_temp - 40.0,
5378 None,
5379 vec![true_density],
5380 )
5381 .unwrap()
5382 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5383 .with_fit_temperature(true)
5384 .with_energy_scale(0.0, 1.0, flight_path);
5385
5386 let error = fit_spectrum_typed(
5387 &InputData::Transmission {
5388 transmission: t_obs,
5389 uncertainty: sigma,
5390 },
5391 &config,
5392 )
5393 .expect_err("transmission plus Poisson/KL must be rejected");
5394 assert!(error.to_string().contains("normalized transmission"));
5395 }
5396
5397 #[test]
5403 fn test_energy_scale_with_temperature_recovers_all_three_counts() {
5404 let data = u238_three_resonances();
5405 let flight_path = 25.0_f64;
5406 let true_density = 0.002;
5407 let true_temp = 340.0;
5408 let true_t0 = 0.6_f64;
5409 let true_l_scale = 1.004_f64;
5410 let flux = 1.0e4_f64;
5411 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5412 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5413 &nominal,
5414 true_t0,
5415 true_l_scale,
5416 flight_path,
5417 )
5418 .unwrap();
5419 let (t_true, _) = synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5420 let open_beam = vec![flux; t_true.len()];
5421 let sample_counts: Vec<f64> = t_true.iter().map(|&t| flux * t).collect();
5422
5423 let config = UnifiedFitConfig::new(
5424 nominal,
5425 vec![data],
5426 vec!["U-238".into()],
5427 true_temp - 40.0,
5428 None,
5429 vec![true_density],
5430 )
5431 .unwrap()
5432 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5433 .with_fit_temperature(true)
5434 .with_energy_scale(0.0, 1.0, flight_path);
5435
5436 let result = fit_spectrum_typed(
5437 &InputData::Counts {
5438 sample_counts,
5439 open_beam_counts: open_beam,
5440 },
5441 &config,
5442 )
5443 .expect("counts joint-Poisson joint fit runs");
5444 assert!(result.converged, "counts joint fit should converge");
5445 let t0 = result.t0_us.expect("t0_us populated");
5446 let ls = result.l_scale.expect("l_scale populated");
5447 let temp = result.temperature_k.expect("temperature_k populated");
5448 assert!(
5449 (t0 - true_t0).abs() < 0.1,
5450 "counts t0: fitted={t0}, true={true_t0}"
5451 );
5452 assert!(
5453 (ls - true_l_scale).abs() / true_l_scale < 2e-3,
5454 "counts l_scale: fitted={ls}, true={true_l_scale}"
5455 );
5456 assert!(
5457 (temp - true_temp).abs() < 5.0,
5458 "counts temperature: fitted={temp}, true={true_temp}"
5459 );
5460 let t_unc = result
5461 .temperature_k_unc
5462 .expect("counts joint fit must report a temperature σ");
5463 assert!(
5464 t_unc.is_finite() && t_unc > 0.0,
5465 "counts joint temperature σ must be finite positive, got {t_unc}"
5466 );
5467 }
5468
5469 #[test]
5487 fn test_energy_scale_temperature_sigma_covariance_plumbing() {
5488 let data = u238_three_resonances();
5489 let flight_path = 25.0_f64;
5490 let true_density = 0.002;
5491 let true_temp = 340.0;
5492 let true_t0 = 0.6_f64;
5493 let true_l_scale = 1.004_f64;
5494 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5495 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5496 &nominal,
5497 true_t0,
5498 true_l_scale,
5499 flight_path,
5500 )
5501 .unwrap();
5502 let (t_clean, sigma) =
5503 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5504 let noise = seeded_gaussian(t_clean.len(), 0x6340_0000_0000_0634);
5507 let t_obs: Vec<f64> = t_clean
5508 .iter()
5509 .zip(sigma.iter())
5510 .zip(noise.iter())
5511 .map(|((&t, &s), &g)| t + s * g)
5512 .collect();
5513
5514 let joint_cfg = UnifiedFitConfig::new(
5518 nominal.clone(),
5519 vec![data.clone()],
5520 vec!["U-238".into()],
5521 true_temp - 40.0,
5522 None,
5523 vec![true_density],
5524 )
5525 .unwrap()
5526 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5527 .with_fit_temperature(true)
5528 .with_energy_scale(0.0, 1.0, flight_path)
5529 .with_tzero_jacobian_method(Some(EnergyScaleJacobianMethod::FiniteDifference));
5530 let joint = fit_spectrum_typed(
5531 &InputData::Transmission {
5532 transmission: t_obs.clone(),
5533 uncertainty: sigma.clone(),
5534 },
5535 &joint_cfg,
5536 )
5537 .unwrap();
5538 let sigma_t_joint = joint.temperature_k_unc.expect("joint σ_T");
5539
5540 let ref_cfg = UnifiedFitConfig::new(
5543 e_true,
5544 vec![data.clone()],
5545 vec!["U-238".into()],
5546 true_temp - 40.0,
5547 None,
5548 vec![true_density],
5549 )
5550 .unwrap()
5551 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5552 .with_fit_temperature(true);
5553 let reference = fit_spectrum_typed(
5554 &InputData::Transmission {
5555 transmission: t_obs,
5556 uncertainty: sigma.clone(),
5557 },
5558 &ref_cfg,
5559 )
5560 .unwrap();
5561 let sigma_t_ref = reference.temperature_k_unc.expect("reference σ_T");
5562
5563 assert!(sigma_t_joint.is_finite() && sigma_t_joint > 0.0);
5571 assert!(sigma_t_ref.is_finite() && sigma_t_ref > 0.0);
5572 assert!(
5573 sigma_t_joint > sigma_t_ref,
5574 "joint σ_T ({sigma_t_joint:.5}) must strictly exceed no-energy-scale \
5575 σ_T ({sigma_t_ref:.5}): the extra t0/L_scale params are not \
5576 orthogonal to temperature"
5577 );
5578
5579 let p = [
5582 joint.densities[0],
5583 joint.temperature_k.unwrap(),
5584 joint.t0_us.unwrap(),
5585 joint.l_scale.unwrap(),
5586 ];
5587 let model = EnergyScaleTransmissionModel::new(
5588 Arc::new(vec![data]),
5589 Arc::new(vec![0]),
5590 Arc::new(vec![1.0]),
5591 true_temp - 40.0,
5592 nominal.clone(),
5593 flight_path,
5594 2,
5595 3,
5596 None,
5597 )
5598 .with_temperature_index(Some(1))
5599 .expect("distinct temperature index")
5600 .with_jacobian_method(EnergyScaleJacobianMethod::FiniteDifference);
5601
5602 let steps = [1e-4 * p[0].max(1e-6), 1e-4 * p[1].max(1.0), 1e-4, 1e-7];
5607 let n_e = nominal.len();
5608 let mut jac = vec![[0.0f64; 4]; n_e];
5609 for (j, &h) in steps.iter().enumerate() {
5610 let mut pp = p;
5611 let mut pm = p;
5612 pp[j] += h;
5613 pm[j] -= h;
5614 let yp = model.evaluate(&pp).unwrap();
5615 let ym = model.evaluate(&pm).unwrap();
5616 for i in 0..n_e {
5617 jac[i][j] = (yp[i] - ym[i]) / (2.0 * h);
5618 }
5619 }
5620 let mut a = vec![vec![0.0f64; 4]; 4];
5622 for i in 0..n_e {
5623 let w = 1.0 / (sigma[i] * sigma[i]);
5624 for r in 0..4 {
5625 for c in 0..4 {
5626 a[r][c] += jac[i][r] * w * jac[i][c];
5627 }
5628 }
5629 }
5630 let cov = invert_dense(&a).expect("covariance invertible");
5631 let sigma_t_recon = (cov[1][1] * joint.reduced_chi_squared).sqrt();
5633 let rel = (sigma_t_recon - sigma_t_joint).abs() / sigma_t_joint;
5634 assert!(
5635 rel < 1e-3,
5636 "independent covariance reconstruction σ_T={sigma_t_recon:.5} must match \
5637 the solver-reported σ_T={sigma_t_joint:.5} to <0.1%, got {rel:.3e} \
5638 (a mismatch means the solver read T's σ from the wrong covariance slot)"
5639 );
5640 }
5641
5642 #[test]
5645 fn test_spectrum_result_corrected_energies() {
5646 let nominal: Vec<f64> = (0..50).map(|i| 5.0 + i as f64 * 0.3).collect();
5647 let flight_path = 25.0;
5648 let (t0, ls) = (0.4_f64, 1.003_f64);
5649 let kl = nereids_physics::resolution::TOF_FACTOR * flight_path;
5655 let expected: Vec<f64> = nominal
5656 .iter()
5657 .map(|&e| {
5658 let tof = kl / e.sqrt();
5659 (kl * ls / (tof - t0)).powi(2)
5660 })
5661 .collect();
5662
5663 let mk = |t0: Option<f64>, ls: Option<f64>, fp: Option<f64>| SpectrumFitResult {
5664 densities: vec![0.001],
5665 uncertainties: None,
5666 reduced_chi_squared: 1.0,
5667 converged: true,
5668 iterations: 1,
5669 temperature_k: None,
5670 temperature_k_unc: None,
5671 anorm: 1.0,
5672 background: [0.0; 3],
5673 back_d: None,
5674 back_f: None,
5675 t0_us: t0,
5676 l_scale: ls,
5677 energy_scale_flight_path_m: fp,
5678 deviance_per_dof: None,
5679 baseline: None,
5680 baseline_e_ref_ev: None,
5681 warnings: Vec::new(),
5682 };
5683
5684 let got = mk(Some(t0), Some(ls), Some(flight_path))
5687 .corrected_energies(&nominal)
5688 .expect("Some when energy scale fitted")
5689 .unwrap();
5690 assert_eq!(got, expected);
5691 assert!(mk(None, None, None).corrected_energies(&nominal).is_none());
5693 assert!(
5694 mk(Some(t0), None, None)
5695 .corrected_energies(&nominal)
5696 .is_none()
5697 );
5698
5699 assert!(
5705 mk(Some(t0), Some(ls), Some(0.0))
5706 .corrected_energies(&nominal)
5707 .unwrap()
5708 .is_err()
5709 );
5710 assert!(
5711 mk(Some(t0), Some(ls), Some(-25.0))
5712 .corrected_energies(&nominal)
5713 .unwrap()
5714 .is_err()
5715 );
5716 assert!(
5717 mk(Some(t0), Some(-1.0), Some(flight_path))
5718 .corrected_energies(&nominal)
5719 .unwrap()
5720 .is_err()
5721 );
5722 assert!(
5723 mk(Some(t0), Some(f64::NAN), Some(flight_path))
5724 .corrected_energies(&nominal)
5725 .unwrap()
5726 .is_err()
5727 );
5728 assert!(
5732 mk(Some(0.0), Some(1.0), Some(flight_path))
5733 .corrected_energies(&[1.0, f64::NAN, 3.0])
5734 .unwrap()
5735 .is_err()
5736 );
5737 assert!(
5738 mk(Some(t0), Some(ls), Some(flight_path))
5739 .corrected_energies(&[0.0, 1.0])
5740 .unwrap()
5741 .is_err()
5742 );
5743 assert!(
5746 mk(Some(t0), Some(ls), Some(flight_path))
5747 .corrected_energies(&[5.0, 4.0, 6.0])
5748 .unwrap()
5749 .is_err()
5750 );
5751 }
5752
5753 #[test]
5755 fn test_energy_scale_returns_fitted_params() {
5756 let data = u238_single_resonance();
5757 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5758 let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);
5759
5760 let config = UnifiedFitConfig::new(
5761 energies,
5762 vec![data],
5763 vec!["U-238".into()],
5764 293.6,
5765 None,
5766 vec![0.001],
5767 )
5768 .unwrap()
5769 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5770 .with_transmission_background(BackgroundConfig::default())
5771 .with_energy_scale(0.0, 1.0, 25.0);
5772
5773 let input = InputData::Transmission {
5774 transmission: t_obs,
5775 uncertainty: sigma,
5776 };
5777 let result = fit_spectrum_typed(&input, &config).unwrap();
5778 assert!(result.converged, "Fit should converge");
5779 assert!(
5780 result.t0_us.is_some(),
5781 "t0_us should be Some when energy-scale is fitted"
5782 );
5783 assert!(
5784 result.l_scale.is_some(),
5785 "l_scale should be Some when energy-scale is fitted"
5786 );
5787 let t0 = result.t0_us.unwrap();
5788 let ls = result.l_scale.unwrap();
5789 assert!(t0.is_finite(), "t0 should be finite, got {t0}");
5791 assert!(ls.is_finite(), "l_scale should be finite, got {ls}");
5792 assert!(t0.abs() < 10.0, "t0 should be within bounds, got {t0}");
5793 assert!(
5794 ls > 0.98 && ls < 1.02,
5795 "l_scale should be within bounds, got {ls}"
5796 );
5797 }
5798
5799 #[test]
5801 fn test_no_energy_scale_returns_none() {
5802 let data = u238_single_resonance();
5803 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5804 let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);
5805
5806 let config = UnifiedFitConfig::new(
5807 energies,
5808 vec![data],
5809 vec!["U-238".into()],
5810 293.6,
5811 None,
5812 vec![0.001],
5813 )
5814 .unwrap()
5815 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5816 .with_transmission_background(BackgroundConfig::default());
5817
5818 let input = InputData::Transmission {
5819 transmission: t_obs,
5820 uncertainty: sigma,
5821 };
5822 let result = fit_spectrum_typed(&input, &config).unwrap();
5823 assert!(
5824 result.t0_us.is_none(),
5825 "t0_us should be None without energy-scale"
5826 );
5827 assert!(
5828 result.l_scale.is_none(),
5829 "l_scale should be None without energy-scale"
5830 );
5831 }
5832
5833 #[test]
5843 fn test_joint_poisson_density_recovery_c_5_98() {
5844 let data = u238_single_resonance();
5845 let true_density = 0.0005_f64;
5846 let c = 5.98_f64;
5847 let lam_ob = 1000.0_f64; let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5849 let (t, _) = synthetic_transmission(&data, true_density, &energies);
5850
5851 let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
5854 let sample_counts: Vec<f64> = t.iter().map(|&ti| c * lam_ob * ti).collect();
5855
5856 let config = UnifiedFitConfig::new(
5857 energies,
5858 vec![data],
5859 vec!["U-238".into()],
5860 0.0,
5861 None,
5862 vec![0.001],
5863 )
5864 .unwrap()
5865 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5866 .with_counts_background(CountsBackgroundConfig {
5867 c,
5868 ..Default::default()
5869 });
5870
5871 let input = InputData::Counts {
5872 sample_counts,
5873 open_beam_counts,
5874 };
5875 let result = fit_spectrum_typed(&input, &config).unwrap();
5876
5877 let d_per_dof = result
5879 .deviance_per_dof
5880 .expect("joint-Poisson solver must populate deviance_per_dof");
5881 assert!(
5882 d_per_dof.is_finite() && d_per_dof >= 0.0,
5883 "deviance_per_dof = {d_per_dof} is not a valid GOF"
5884 );
5885 assert!(
5889 d_per_dof < 0.5,
5890 "noise-free expected-counts fit should give D/dof ≈ 0, got {d_per_dof}"
5891 );
5892 let rel_bias = (result.densities[0] - true_density) / true_density;
5894 assert!(
5895 rel_bias.abs() < 0.05,
5896 "density bias {rel_bias} > 5%: fitted={} truth={true_density}",
5897 result.densities[0]
5898 );
5899 assert!((result.reduced_chi_squared - d_per_dof).abs() < 1e-12);
5901 }
5902
5903 #[test]
5907 fn test_joint_poisson_rejects_alpha_fit() {
5908 let data = u238_single_resonance();
5909 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
5910 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
5911 let open_beam_counts: Vec<f64> = vec![500.0; energies.len()];
5912 let sample_counts: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
5913
5914 let config = UnifiedFitConfig::new(
5915 energies,
5916 vec![data],
5917 vec!["U-238".into()],
5918 0.0,
5919 None,
5920 vec![0.001],
5921 )
5922 .unwrap()
5923 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5924 .with_counts_background(CountsBackgroundConfig {
5925 fit_alpha_1: true,
5926 c: 1.0,
5927 ..Default::default()
5928 });
5929
5930 let input = InputData::Counts {
5931 sample_counts,
5932 open_beam_counts,
5933 };
5934 let err = fit_spectrum_typed(&input, &config).unwrap_err();
5935 let msg = err.to_string();
5936 assert!(
5937 msg.contains("fit_alpha_1") || msg.contains("alpha_1"),
5938 "expected alpha_1 rejection message, got: {msg}"
5939 );
5940 }
5941
5942 #[test]
5945 fn test_transmission_poisson_kl_rejects_fractional_data() {
5946 let data = u238_single_resonance();
5947 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
5948 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
5949
5950 let config = UnifiedFitConfig::new(
5951 energies,
5952 vec![data],
5953 vec!["U-238".into()],
5954 0.0,
5955 None,
5956 vec![0.001],
5957 )
5958 .unwrap()
5959 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
5960
5961 let input = InputData::Transmission {
5962 transmission: t,
5963 uncertainty: u,
5964 };
5965 let error = fit_spectrum_typed(&input, &config)
5966 .expect_err("transmission plus Poisson/KL must be rejected");
5967 assert!(error.to_string().contains("Poisson"));
5968 }
5969
5970 #[test]
5974 fn counts_resolution_requires_exact_count_response() {
5975 use nereids_physics::resolution::{
5976 ResolutionFunction, ResolutionParams, TabulatedResolution,
5977 };
5978
5979 let data = u238_single_resonance();
5980 let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
5981 let (sample_counts, open_beam_counts) = synthetic_counts(&data, 0.0005, &energies, 1000.0);
5982 let tab_text = "header\n---\n\
5983 5.0 0.0\n\
5984 -0.01 0.0\n\
5985 0.0 1.0\n\
5986 0.01 0.0\n\
5987 \n\
5988 20.0 0.0\n\
5989 -0.02 0.0\n\
5990 0.0 1.0\n\
5991 0.02 0.0\n";
5992 let resolutions = [
5993 (
5994 "gaussian",
5995 ResolutionFunction::Gaussian(ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap()),
5996 ),
5997 (
5998 "tabulated",
5999 ResolutionFunction::Tabulated(Arc::new(
6000 TabulatedResolution::from_text(tab_text, 25.0).unwrap(),
6001 )),
6002 ),
6003 ];
6004 let solvers = [
6005 ("auto", SolverConfig::Auto),
6006 ("kl", SolverConfig::PoissonKL(PoissonConfig::default())),
6007 ("lm", SolverConfig::LevenbergMarquardt(LmConfig::default())),
6008 ];
6009
6010 for (resolution_name, resolution) in resolutions {
6011 for (solver_name, solver) in &solvers {
6012 let config = UnifiedFitConfig::new(
6013 energies.clone(),
6014 vec![data.clone()],
6015 vec!["U-238".into()],
6016 293.6,
6017 Some(resolution.clone()),
6018 vec![0.001],
6019 )
6020 .unwrap()
6021 .with_solver(solver.clone());
6022 let inputs = [
6023 (
6024 "counts",
6025 InputData::Counts {
6026 sample_counts: sample_counts.clone(),
6027 open_beam_counts: open_beam_counts.clone(),
6028 },
6029 ),
6030 (
6031 "counts-with-nuisance",
6032 InputData::CountsWithNuisance {
6033 sample_counts: sample_counts.clone(),
6034 flux: open_beam_counts.clone(),
6035 background: vec![0.0; energies.len()],
6036 },
6037 ),
6038 ];
6039
6040 for (input_name, input) in inputs {
6041 let err = match fit_spectrum_typed(&input, &config) {
6042 Ok(_) => panic!(
6043 "{input_name} + {solver_name} + {resolution_name} resolution \
6044 returned a scientifically invalid fit"
6045 ),
6046 Err(err) => err,
6047 };
6048 let msg = err.to_string();
6049 assert!(
6050 msg.contains("counts")
6051 && msg.contains("instrument resolution")
6052 && msg.contains("separate-arm model"),
6053 "{input_name} + {solver_name} + {resolution_name}: expected \
6054 physical counts-response rejection, got: {msg}"
6055 );
6056 }
6057 }
6058 }
6059 }
6060
6061 fn run_exact_closed_loop(timing_offset_us: f64) {
6069 use nereids_physics::counts_response::two_arm_count_response;
6070 use nereids_physics::resolution::{ResolutionFunction, TabulatedResolution};
6071
6072 let data = u238_single_resonance();
6073 let true_density = 5.0e-4;
6074 let energies: Vec<f64> = (0..80).map(|i| 5.0 + i as f64 * 3.0 / 79.0).collect();
6075 let (true_transmission, _) = synthetic_transmission(&data, true_density, &energies);
6076 let response = ResolutionFunction::Tabulated(Arc::new(
6077 TabulatedResolution::from_kernels(
6078 vec![6.5],
6079 vec![(vec![-3.0, 0.0, 3.0], vec![0.0, 1.0, 0.0])],
6080 25.0,
6081 )
6082 .expect("valid detector-time response"),
6083 ));
6084 let source: Vec<f64> = energies
6085 .iter()
6086 .enumerate()
6087 .map(|(index, _)| 4.0e4 * (1.0 + 0.4 * index as f64 / 79.0))
6088 .collect();
6089 let detector_edges: Vec<f64> = (0..102)
6092 .map(|i| 630.0 + timing_offset_us + i as f64 * 2.5)
6093 .collect();
6094 let expected = two_arm_count_response(
6095 &energies,
6096 &source,
6097 &true_transmission,
6098 &detector_edges,
6099 timing_offset_us,
6100 &response,
6101 )
6102 .expect("valid synthetic count response");
6103 let fixed_anorm = 0.98;
6104 let fixed_back_a = 0.01;
6105 let fixed_back_b = 0.002;
6106 let sample_counts: Vec<f64> = expected
6107 .open_beam
6108 .iter()
6109 .zip(&expected.sample)
6110 .zip(detector_edges.windows(2))
6111 .map(|((&open, &sample), edges)| {
6112 if open == 0.0 {
6113 0.0
6114 } else {
6115 let measured_energy =
6125 tof_to_energy(0.5 * (edges[0] + edges[1]) - timing_offset_us, 25.0);
6126 open * (fixed_anorm * sample / open
6127 + fixed_back_a
6128 + fixed_back_b / measured_energy.sqrt())
6129 }
6130 })
6131 .collect();
6132
6133 let config = UnifiedFitConfig::new(
6134 energies,
6135 vec![data],
6136 vec!["U-238".into()],
6137 0.0,
6138 Some(response),
6139 vec![2.0e-4],
6140 )
6141 .unwrap()
6142 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6143 max_iter: 400,
6144 ..Default::default()
6145 }))
6146 .with_exact_count_response(ExactCountResponseConfig {
6147 incident_fluence_weights: source,
6148 detector_time_edges_us: detector_edges,
6149 timing_offset_us,
6150 })
6151 .with_transmission_background(BackgroundConfig {
6152 anorm_init: fixed_anorm,
6153 back_a_init: fixed_back_a,
6154 back_b_init: fixed_back_b,
6155 back_c_init: 0.0,
6156 back_d_init: 0.01,
6160 back_f_init: 1.0,
6161 fit_anorm: false,
6162 fit_back_a: false,
6163 fit_back_b: false,
6164 fit_back_c: false,
6165 fit_back_d: false,
6166 fit_back_f: false,
6167 });
6168 let result = fit_spectrum_typed(
6169 &InputData::Counts {
6170 sample_counts,
6171 open_beam_counts: expected.open_beam,
6172 },
6173 &config,
6174 )
6175 .expect("exact resolved count fit");
6176
6177 assert!(result.converged, "fit did not converge");
6178 assert!(
6179 (result.densities[0] - true_density).abs() < 2.0e-6,
6180 "recovered density {} differs from truth {true_density}",
6181 result.densities[0]
6182 );
6183 assert!(result.deviance_per_dof.unwrap_or(f64::INFINITY) < 1.0e-8);
6184 }
6185
6186 #[test]
6187 fn exact_resolved_counts_closed_loop_recovers_density() {
6188 run_exact_closed_loop(0.0);
6189 }
6190
6191 #[test]
6195 fn exact_resolved_counts_closed_loop_with_timing_offset() {
6196 run_exact_closed_loop(5.0);
6197 }
6198
6199 #[test]
6205 fn exact_resolved_counts_recovers_free_normalization() {
6206 use nereids_physics::counts_response::two_arm_count_response;
6207 use nereids_physics::resolution::{ResolutionFunction, TabulatedResolution};
6208
6209 let data = u238_single_resonance();
6210 let true_density = 5.0e-4;
6211 let true_anorm = 0.97;
6212 let true_back_a = 0.015;
6213 let energies: Vec<f64> = (0..80).map(|i| 5.0 + i as f64 * 3.0 / 79.0).collect();
6214 let (true_transmission, _) = synthetic_transmission(&data, true_density, &energies);
6215 let response = ResolutionFunction::Tabulated(Arc::new(
6216 TabulatedResolution::from_kernels(
6217 vec![6.5],
6218 vec![(vec![-3.0, 0.0, 3.0], vec![0.0, 1.0, 0.0])],
6219 25.0,
6220 )
6221 .expect("valid detector-time response"),
6222 ));
6223 let source: Vec<f64> = energies
6224 .iter()
6225 .enumerate()
6226 .map(|(index, _)| 4.0e4 * (1.0 + 0.4 * index as f64 / 79.0))
6227 .collect();
6228 let detector_edges: Vec<f64> = (0..102).map(|i| 630.0 + i as f64 * 2.5).collect();
6229 let expected = two_arm_count_response(
6230 &energies,
6231 &source,
6232 &true_transmission,
6233 &detector_edges,
6234 0.0,
6235 &response,
6236 )
6237 .expect("valid synthetic count response");
6238 let sample_counts: Vec<f64> = expected
6239 .open_beam
6240 .iter()
6241 .zip(&expected.sample)
6242 .map(|(&open, &sample)| {
6243 if open == 0.0 {
6244 0.0
6245 } else {
6246 open * (true_anorm * sample / open + true_back_a)
6247 }
6248 })
6249 .collect();
6250
6251 let config = UnifiedFitConfig::new(
6252 energies,
6253 vec![data],
6254 vec!["U-238".into()],
6255 0.0,
6256 Some(response),
6257 vec![2.0e-4],
6258 )
6259 .unwrap()
6260 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6261 max_iter: 600,
6262 ..Default::default()
6263 }))
6264 .with_exact_count_response(ExactCountResponseConfig {
6265 incident_fluence_weights: source,
6266 detector_time_edges_us: detector_edges,
6267 timing_offset_us: 0.0,
6268 })
6269 .with_transmission_background(BackgroundConfig {
6270 anorm_init: 1.0, back_a_init: 0.0, back_b_init: 0.0,
6273 back_c_init: 0.0,
6274 back_d_init: 0.01,
6275 back_f_init: 1.0,
6276 fit_anorm: true,
6277 fit_back_a: true,
6278 fit_back_b: false,
6279 fit_back_c: false,
6280 fit_back_d: false,
6281 fit_back_f: false,
6282 });
6283 let result = fit_spectrum_typed(
6284 &InputData::Counts {
6285 sample_counts,
6286 open_beam_counts: expected.open_beam,
6287 },
6288 &config,
6289 )
6290 .expect("exact fit with free normalization");
6291
6292 assert!(result.converged, "fit did not converge");
6293 assert!(
6294 (result.densities[0] - true_density).abs() / true_density < 1.0e-3,
6295 "density: fitted={}, true={true_density}",
6296 result.densities[0]
6297 );
6298 assert!(
6299 (result.anorm - true_anorm).abs() < 1.0e-3,
6300 "anorm: fitted={}, true={true_anorm}",
6301 result.anorm
6302 );
6303 assert!(
6304 (result.background[0] - true_back_a).abs() < 1.0e-3,
6305 "back_a: fitted={}, true={true_back_a}",
6306 result.background[0]
6307 );
6308 assert!(result.deviance_per_dof.unwrap_or(f64::INFINITY) < 1.0e-6);
6309 }
6310
6311 #[test]
6318 fn exact_resolved_counts_recovers_free_baseline_tilt() {
6319 use nereids_physics::counts_response::two_arm_count_response;
6320 use nereids_physics::resolution::{ResolutionFunction, TabulatedResolution};
6321
6322 let data = u238_single_resonance();
6323 let true_density = 5.0e-4;
6324 let true_b1 = 0.02;
6325 let energies: Vec<f64> = (0..80).map(|i| 5.0 + i as f64 * 3.0 / 79.0).collect();
6326 let (true_transmission, _) = synthetic_transmission(&data, true_density, &energies);
6327 let e_ref =
6331 nereids_fitting::transmission_model::baseline_reference_energy_active(&energies, None);
6332 let baselined: Vec<f64> = energies
6333 .iter()
6334 .zip(&true_transmission)
6335 .map(|(&e, &t)| (1.0 + true_b1 * (e / e_ref).ln()) * t)
6336 .collect();
6337 let response = ResolutionFunction::Tabulated(Arc::new(
6338 TabulatedResolution::from_kernels(
6339 vec![6.5],
6340 vec![(vec![-3.0, 0.0, 3.0], vec![0.0, 1.0, 0.0])],
6341 25.0,
6342 )
6343 .expect("valid detector-time response"),
6344 ));
6345 let source: Vec<f64> = energies
6346 .iter()
6347 .enumerate()
6348 .map(|(index, _)| 4.0e4 * (1.0 + 0.4 * index as f64 / 79.0))
6349 .collect();
6350 let detector_edges: Vec<f64> = (0..102).map(|i| 630.0 + i as f64 * 2.5).collect();
6351 let expected = two_arm_count_response(
6352 &energies,
6353 &source,
6354 &baselined,
6355 &detector_edges,
6356 0.0,
6357 &response,
6358 )
6359 .expect("valid synthetic count response");
6360
6361 let config = UnifiedFitConfig::new(
6362 energies,
6363 vec![data],
6364 vec!["U-238".into()],
6365 0.0,
6366 Some(response),
6367 vec![2.0e-4],
6368 )
6369 .unwrap()
6370 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6371 max_iter: 600,
6372 ..Default::default()
6373 }))
6374 .with_exact_count_response(ExactCountResponseConfig {
6375 incident_fluence_weights: source,
6376 detector_time_edges_us: detector_edges,
6377 timing_offset_us: 0.0,
6378 })
6379 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
6380 b0_init: 1.0,
6381 b1_init: 0.0, b2_init: 0.0,
6383 fit_b0: false,
6384 fit_b1: true,
6385 fit_b2: false,
6386 ..Default::default()
6387 });
6388 let result = fit_spectrum_typed(
6389 &InputData::Counts {
6390 sample_counts: expected.sample,
6391 open_beam_counts: expected.open_beam,
6392 },
6393 &config,
6394 )
6395 .expect("exact fit with free baseline tilt");
6396
6397 assert!(result.converged, "fit did not converge");
6398 assert!(
6399 (result.densities[0] - true_density).abs() / true_density < 1.0e-3,
6400 "density: fitted={}, true={true_density}",
6401 result.densities[0]
6402 );
6403 let baseline = result.baseline.expect("baseline coefficients reported");
6404 assert!(
6405 (baseline[1] - true_b1).abs() < 1.0e-3,
6406 "b1: fitted={}, true={true_b1}",
6407 baseline[1]
6408 );
6409 assert!(result.deviance_per_dof.unwrap_or(f64::INFINITY) < 1.0e-6);
6410 }
6411
6412 #[test]
6418 fn exact_route_validation_arms_reject_each_misuse() {
6419 use nereids_physics::resolution::{ResolutionFunction, TOF_FACTOR, TabulatedResolution};
6420
6421 let make_response = || {
6422 ResolutionFunction::Tabulated(Arc::new(
6423 TabulatedResolution::from_kernels(
6424 vec![25.0],
6425 vec![(vec![-1.0, 0.0, 1.0], vec![0.0, 1.0, 0.0])],
6426 25.0,
6427 )
6428 .expect("valid triangle response"),
6429 ))
6430 };
6431 let arrival = TOF_FACTOR * 25.0 / 25.0_f64.sqrt();
6432 let data = u238_single_resonance();
6433 let base_config = |resolution: Option<ResolutionFunction>| {
6434 UnifiedFitConfig::new(
6435 vec![25.0],
6436 vec![data.clone()],
6437 vec!["U-238".into()],
6438 0.0,
6439 resolution,
6440 vec![1.0e-4],
6441 )
6442 .unwrap()
6443 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6444 };
6445 let exact = |edges: Vec<f64>| ExactCountResponseConfig {
6446 incident_fluence_weights: vec![100.0],
6447 detector_time_edges_us: edges,
6448 timing_offset_us: 0.0,
6449 };
6450 let two_bin_edges = vec![arrival - 1.0, arrival, arrival + 1.0];
6451 let counts_2 = InputData::Counts {
6452 sample_counts: vec![10.0, 10.0],
6453 open_beam_counts: vec![50.0, 50.0],
6454 };
6455
6456 let err = fit_spectrum_typed(
6459 &counts_2,
6460 &base_config(Some(make_response())).with_exact_count_response(exact(vec![
6461 arrival + 10.0,
6462 arrival + 11.0,
6463 arrival + 12.0,
6464 ])),
6465 )
6466 .expect_err("occupied dead bin must be a hard error");
6467 assert!(
6468 err.to_string()
6469 .contains("zero detector response in occupied detector bin"),
6470 "{err}"
6471 );
6472
6473 let err = fit_spectrum_typed(
6480 &InputData::CountsWithNuisance {
6481 sample_counts: vec![10.0, 10.0],
6482 flux: vec![50.0, 50.0],
6483 background: vec![5.0, 5.0],
6484 },
6485 &base_config(Some(make_response())).with_exact_count_response(exact(vec![
6486 arrival + 10.0,
6487 arrival + 11.0,
6488 arrival + 12.0,
6489 ])),
6490 )
6491 .expect_err("dropping every bin must be reported, not fitted");
6492 let message = err.to_string();
6493 assert!(
6494 message.contains("hold only background and cannot be fitted"),
6495 "a declared background must reclassify the dead bin rather than \
6496 rejecting the source: {message}"
6497 );
6498 assert!(
6499 !message.contains("zero detector response in occupied detector bin"),
6500 "the no-background rejection must not fire once a background is \
6501 declared: {message}"
6502 );
6503
6504 let err = fit_spectrum_typed(
6509 &InputData::CountsWithNuisance {
6510 sample_counts: vec![10.0, 10.0],
6511 flux: vec![50.0, 50.0],
6512 background: vec![0.0, 5.0],
6513 },
6514 &base_config(Some(make_response())).with_exact_count_response(exact(vec![
6515 arrival + 10.0,
6516 arrival + 11.0,
6517 arrival + 12.0,
6518 ])),
6519 )
6520 .expect_err("a dead bin with no background of its own is still unexplained");
6521 assert!(
6522 err.to_string()
6523 .contains("zero detector response in occupied detector bin 0"),
6524 "the bin without a local background must be the one reported: {err}"
6525 );
6526
6527 let err = fit_spectrum_typed(
6529 &InputData::Transmission {
6530 transmission: vec![0.9],
6531 uncertainty: vec![0.01],
6532 },
6533 &base_config(Some(make_response()))
6534 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6535 .with_exact_count_response(exact(two_bin_edges.clone())),
6536 )
6537 .expect_err("exact config on transmission must be rejected");
6538 assert!(
6539 err.to_string()
6540 .contains("cannot be attached to normalized transmission"),
6541 "{err}"
6542 );
6543
6544 let err = fit_spectrum_typed(
6546 &counts_2,
6547 &base_config(None).with_exact_count_response(exact(two_bin_edges.clone())),
6548 )
6549 .expect_err("exact config without resolution must be rejected");
6550 assert!(
6551 err.to_string()
6552 .contains("requires an instrument resolution model"),
6553 "{err}"
6554 );
6555
6556 let err = fit_spectrum_typed(
6558 &counts_2,
6559 &base_config(Some(make_response())).with_exact_count_response(
6560 ExactCountResponseConfig {
6561 incident_fluence_weights: vec![100.0, 100.0],
6562 detector_time_edges_us: two_bin_edges.clone(),
6563 timing_offset_us: 0.0,
6564 },
6565 ),
6566 )
6567 .expect_err("fluence length mismatch must be rejected");
6568 assert!(
6569 err.to_string().contains("incident_fluence_weights length"),
6570 "{err}"
6571 );
6572
6573 let err = fit_spectrum_typed(
6575 &counts_2,
6576 &base_config(Some(make_response()))
6577 .with_exact_count_response(exact(vec![arrival - 1.0, arrival])),
6578 )
6579 .expect_err("edge/bin length mismatch must be rejected");
6580 assert!(
6581 err.to_string().contains("detector_time_edges_us length"),
6582 "{err}"
6583 );
6584
6585 let err = fit_spectrum_typed(
6587 &counts_2,
6588 &base_config(Some(make_response())).with_exact_count_response(
6589 ExactCountResponseConfig {
6590 incident_fluence_weights: vec![100.0],
6591 detector_time_edges_us: two_bin_edges.clone(),
6592 timing_offset_us: f64::NAN,
6593 },
6594 ),
6595 )
6596 .expect_err("NaN timing offset must be rejected");
6597 assert!(
6598 err.to_string().contains("timing_offset_us must be finite"),
6599 "{err}"
6600 );
6601
6602 let err = fit_spectrum_typed(
6604 &counts_2,
6605 &base_config(Some(make_response()))
6606 .with_exact_count_response(exact(two_bin_edges.clone()))
6607 .with_energy_scale(0.0, 1.0, 25.0),
6608 )
6609 .expect_err("energy-scale through the exact response must be rejected");
6610 assert!(
6611 err.to_string()
6612 .contains("not yet connected to the exact two-arm"),
6613 "{err}"
6614 );
6615
6616 let err = fit_spectrum_typed(
6618 &counts_2,
6619 &base_config(Some(make_response()))
6620 .with_exact_count_response(exact(two_bin_edges.clone()))
6621 .with_fit_energy_range(Some((10.0, 30.0)))
6622 .unwrap(),
6623 )
6624 .expect_err("fit_energy_range through the exact response must be rejected");
6625 assert!(
6626 err.to_string()
6627 .contains("fit_energy_range is not yet supported"),
6628 "{err}"
6629 );
6630
6631 let sample = ndarray::Array3::from_elem((1, 2, 2), 10.0);
6633 let open_beam = ndarray::Array3::from_elem((1, 2, 2), 50.0);
6634 let err = crate::spatial::spatial_map_typed(
6635 &crate::spatial::InputData3D::Counts {
6636 sample_counts: sample.view(),
6637 open_beam_counts: open_beam.view(),
6638 },
6639 &base_config(Some(make_response())).with_exact_count_response(exact(two_bin_edges)),
6640 None,
6641 None,
6642 None,
6643 )
6644 .expect_err("spatial exact-count mapping must be rejected");
6645 assert!(
6646 err.to_string()
6647 .contains("single-spectrum count fitter only"),
6648 "{err}"
6649 );
6650 }
6651
6652 fn run_hand_computed_anchor(timing_offset_us: f64) {
6668 use nereids_physics::resolution::{ResolutionFunction, TOF_FACTOR, TabulatedResolution};
6669
6670 let flight_path_m = 25.0_f64;
6671 let arrival_0 = TOF_FACTOR * flight_path_m / 25.0_f64.sqrt();
6672 let energy_1 = (TOF_FACTOR * flight_path_m / (arrival_0 + 1.0)).powi(2);
6673 let energies = vec![25.0, energy_1];
6674 let response = ResolutionFunction::Tabulated(Arc::new(
6675 TabulatedResolution::from_kernels(
6676 vec![25.0],
6677 vec![(vec![-1.0, 0.0, 1.0], vec![0.0, 1.0, 0.0])],
6678 flight_path_m,
6679 )
6680 .expect("valid triangle response"),
6681 ));
6682 let base = timing_offset_us + arrival_0;
6689 let detector_edges = vec![base - 1.0, base, base + 1.0, base + 2.0];
6690 let open_beam_counts = vec![50.0, 150.0, 100.0];
6691 let t_eff = [0.2_f64, 0.6, 0.8];
6692
6693 let fixed_anorm = 1.0_f64;
6696 let fixed_back_b = 0.05_f64;
6697 let sample_counts: Vec<f64> = detector_edges
6698 .windows(2)
6699 .zip(t_eff)
6700 .zip(&open_beam_counts)
6701 .map(|((edges, t), &open)| {
6702 let corrected_tof_us = 0.5 * (edges[0] + edges[1]) - timing_offset_us;
6703 let pseudo_energy = (TOF_FACTOR * flight_path_m / corrected_tof_us).powi(2);
6704 open * (fixed_anorm * t + fixed_back_b / pseudo_energy.sqrt())
6705 })
6706 .collect();
6707
6708 let true_density = 1.0e-3;
6710 let data = u238_single_resonance();
6711 let xs: Vec<Vec<f64>> = vec![vec![
6712 -(0.2_f64.ln()) / true_density,
6713 -(0.8_f64.ln()) / true_density,
6714 ]];
6715
6716 let config = UnifiedFitConfig::new(
6717 energies,
6718 vec![data],
6719 vec!["U-238".into()],
6720 0.0,
6721 Some(response),
6722 vec![5.0e-4], )
6724 .unwrap();
6725 let xs = table_on_data_grid(&config, xs);
6726 let config = config
6727 .with_precomputed_cross_sections(xs)
6728 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6729 max_iter: 200,
6730 ..Default::default()
6731 }))
6732 .with_exact_count_response(ExactCountResponseConfig {
6733 incident_fluence_weights: vec![100.0, 200.0],
6734 detector_time_edges_us: detector_edges,
6735 timing_offset_us,
6736 })
6737 .with_transmission_background(BackgroundConfig {
6738 anorm_init: fixed_anorm,
6739 back_a_init: 0.0,
6740 back_b_init: fixed_back_b,
6741 back_c_init: 0.0,
6742 back_d_init: 0.0,
6743 back_f_init: 1.0,
6744 fit_anorm: false,
6745 fit_back_a: false,
6746 fit_back_b: false,
6747 fit_back_c: false,
6748 fit_back_d: false,
6749 fit_back_f: false,
6750 });
6751 let result = fit_spectrum_typed(
6752 &InputData::Counts {
6753 sample_counts,
6754 open_beam_counts,
6755 },
6756 &config,
6757 )
6758 .expect("hand-anchored exact fit");
6759
6760 assert!(result.converged, "fit did not converge");
6761 assert!(
6762 (result.densities[0] - true_density).abs() / true_density < 1.0e-6,
6763 "recovered density {} differs from hand-derived truth {true_density}",
6764 result.densities[0]
6765 );
6766 assert!(result.deviance_per_dof.unwrap_or(f64::INFINITY) < 1.0e-10);
6767 }
6768
6769 #[test]
6770 fn exact_route_matches_hand_computed_detector_counts() {
6771 run_hand_computed_anchor(0.0);
6772 }
6773
6774 #[test]
6784 fn exact_route_tolerates_empty_pre_trigger_bins_with_background() {
6785 use nereids_physics::counts_response::two_arm_count_response;
6786 use nereids_physics::resolution::{ResolutionFunction, TabulatedResolution};
6787
6788 let timing_offset_us = 5.0_f64;
6789 let data = u238_single_resonance();
6790 let true_density = 5.0e-4;
6791 let energies: Vec<f64> = (0..60).map(|i| 5.0 + i as f64 * 3.0 / 59.0).collect();
6792 let (true_transmission, _) = synthetic_transmission(&data, true_density, &energies);
6793 let response = ResolutionFunction::Tabulated(Arc::new(
6794 TabulatedResolution::from_kernels(
6795 vec![6.5],
6796 vec![(vec![-3.0, 0.0, 3.0], vec![0.0, 1.0, 0.0])],
6797 25.0,
6798 )
6799 .expect("valid detector-time response"),
6800 ));
6801 let source: Vec<f64> = vec![4.0e4; energies.len()];
6802 let detector_edges: Vec<f64> = (0..340).map(|i| i as f64 * 2.5).collect();
6805 let expected = two_arm_count_response(
6806 &energies,
6807 &source,
6808 &true_transmission,
6809 &detector_edges,
6810 timing_offset_us,
6811 &response,
6812 )
6813 .expect("valid synthetic count response");
6814 assert!(
6815 expected.open_beam[0] == 0.0 && expected.sample[0] == 0.0,
6816 "fixture must have an empty pre-trigger bin 0"
6817 );
6818
6819 let config = UnifiedFitConfig::new(
6820 energies,
6821 vec![data],
6822 vec!["U-238".into()],
6823 0.0,
6824 Some(response),
6825 vec![2.0e-4],
6826 )
6827 .unwrap()
6828 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6829 max_iter: 400,
6830 ..Default::default()
6831 }))
6832 .with_exact_count_response(ExactCountResponseConfig {
6833 incident_fluence_weights: source,
6834 detector_time_edges_us: detector_edges,
6835 timing_offset_us,
6836 })
6837 .with_transmission_background(BackgroundConfig {
6839 anorm_init: 1.0,
6840 back_a_init: 0.0,
6841 back_b_init: 0.0,
6842 back_c_init: 0.0,
6843 back_d_init: 0.0,
6844 back_f_init: 1.0,
6845 fit_anorm: false,
6846 fit_back_a: true,
6847 fit_back_b: false,
6848 fit_back_c: false,
6849 fit_back_d: false,
6850 fit_back_f: false,
6851 });
6852 let result = fit_spectrum_typed(
6853 &InputData::Counts {
6854 sample_counts: expected.sample,
6855 open_beam_counts: expected.open_beam,
6856 },
6857 &config,
6858 )
6859 .expect("empty pre-trigger bins must not reject a background fit");
6860
6861 assert!(result.converged, "fit did not converge");
6862 assert!(
6863 (result.densities[0] - true_density).abs() / true_density < 1.0e-3,
6864 "density: fitted={}, true={true_density}",
6865 result.densities[0]
6866 );
6867 }
6868
6869 #[test]
6880 fn exact_resolved_counts_closed_loop_with_ikeda_carpenter() {
6881 use nereids_physics::counts_response::two_arm_count_response;
6882 use nereids_physics::ikeda_carpenter::{
6883 IkedaCarpenter, IkedaCarpenterParams, SynthesisGrid,
6884 };
6885 use nereids_physics::resolution::{ResolutionFunction, TOF_FACTOR};
6886
6887 let flight_path_m = 25.0_f64;
6888 let data = u238_single_resonance();
6889 let true_density = 5.0e-4;
6890 let energies: Vec<f64> = (0..60).map(|i| 5.0 + i as f64 * 3.0 / 59.0).collect();
6891 let (true_transmission, _) = synthetic_transmission(&data, true_density, &energies);
6892 let response = ResolutionFunction::IkedaCarpenter(Arc::new(
6893 IkedaCarpenter::new(
6894 IkedaCarpenterParams::constant(2.0, 0.1, 0.0),
6895 flight_path_m,
6896 &SynthesisGrid::new(4.0, 9.0),
6897 )
6898 .expect("valid prompt-only IC model"),
6899 ));
6900 let source: Vec<f64> = vec![4.0e4; energies.len()];
6901 let first_arrival = TOF_FACTOR * flight_path_m / 8.0_f64.sqrt();
6904 let detector_edges: Vec<f64> = (0..160).map(|i| first_arrival + i as f64 * 2.0).collect();
6905 let expected = two_arm_count_response(
6906 &energies,
6907 &source,
6908 &true_transmission,
6909 &detector_edges,
6910 0.0,
6911 &response,
6912 )
6913 .expect("valid IC synthetic count response");
6914 assert!(
6915 expected.open_beam.iter().any(|&v| v > 0.0),
6916 "IC fixture must deposit counts in the window"
6917 );
6918
6919 let config = UnifiedFitConfig::new(
6920 energies,
6921 vec![data],
6922 vec!["U-238".into()],
6923 0.0,
6924 Some(response),
6925 vec![2.0e-4],
6926 )
6927 .unwrap()
6928 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
6929 max_iter: 400,
6930 ..Default::default()
6931 }))
6932 .with_exact_count_response(ExactCountResponseConfig {
6933 incident_fluence_weights: source,
6934 detector_time_edges_us: detector_edges,
6935 timing_offset_us: 0.0,
6936 });
6937 let result = fit_spectrum_typed(
6938 &InputData::Counts {
6939 sample_counts: expected.sample,
6940 open_beam_counts: expected.open_beam,
6941 },
6942 &config,
6943 )
6944 .expect("exact resolved count fit with an IC response");
6945
6946 assert!(result.converged, "IC exact fit did not converge");
6947 assert!(
6948 (result.densities[0] - true_density).abs() / true_density < 1.0e-3,
6949 "density: fitted={}, true={true_density}",
6950 result.densities[0]
6951 );
6952 }
6953
6954 #[test]
6965 fn exact_route_rejects_occupied_pre_trigger_bin_with_background() {
6966 use nereids_physics::counts_response::DetectorBinResponseMatrix;
6967 use nereids_physics::resolution::{ResolutionFunction, TOF_FACTOR, TabulatedResolution};
6968
6969 let flight_path_m = 25.0_f64;
6970 let true_energy = 32_670.0_f64;
6971 let tof_us = TOF_FACTOR * flight_path_m / true_energy.sqrt();
6972 let timing_offset_us = 30.0_f64;
6973 let arrival = timing_offset_us + tof_us; assert!(
6975 tof_us < 20.0,
6976 "fixture needs a flight time inside the kernel half-width, got {tof_us}"
6977 );
6978 let response = ResolutionFunction::Tabulated(Arc::new(
6979 TabulatedResolution::from_kernels(
6980 vec![true_energy],
6981 vec![(vec![-20.0, 0.0, 20.0], vec![0.0, 1.0, 0.0])],
6982 flight_path_m,
6983 )
6984 .expect("valid wide triangle response"),
6985 ));
6986 let detector_edges = vec![
6989 arrival - 22.0,
6990 arrival - 14.0,
6991 arrival - 6.0,
6992 arrival + 2.0,
6993 arrival + 10.0,
6994 ];
6995 let matrix = DetectorBinResponseMatrix::new(
6998 &[true_energy],
6999 &detector_edges,
7000 timing_offset_us,
7001 &response,
7002 )
7003 .expect("valid response matrix");
7004 assert!(
7005 matrix.probability(0, 0) > 0.0,
7006 "fixture bin 0 must have nonzero detector response"
7007 );
7008 assert!(
7009 0.5 * (detector_edges[0] + detector_edges[1]) < timing_offset_us,
7010 "fixture bin 0 must precede the timing offset"
7011 );
7012
7013 let config = UnifiedFitConfig::new(
7014 vec![true_energy],
7015 vec![u238_single_resonance()],
7016 vec!["U-238".into()],
7017 0.0,
7018 Some(response),
7019 vec![1.0e-4],
7020 )
7021 .unwrap()
7022 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7023 .with_exact_count_response(ExactCountResponseConfig {
7024 incident_fluence_weights: vec![100.0],
7025 detector_time_edges_us: detector_edges,
7026 timing_offset_us,
7027 })
7028 .with_transmission_background(BackgroundConfig {
7029 fit_anorm: false,
7030 fit_back_a: true,
7031 ..BackgroundConfig::default()
7032 });
7033 let err = fit_spectrum_typed(
7034 &InputData::Counts {
7035 sample_counts: vec![5.0, 20.0, 20.0, 5.0],
7036 open_beam_counts: vec![10.0, 50.0, 50.0, 10.0],
7037 },
7038 &config,
7039 )
7040 .expect_err("occupied bin without a physical energy must fail closed");
7041 let msg = err.to_string();
7042 assert!(
7043 msg.contains("carries observed counts"),
7044 "expected the pre-trigger occupied-bin rejection, got: {msg}"
7045 );
7046 }
7047
7048 #[test]
7053 fn exact_route_matches_hand_computed_counts_with_timing_offset() {
7054 run_hand_computed_anchor(7.5);
7055 }
7056
7057 #[test]
7067 fn test_joint_poisson_with_transmission_background() {
7068 let data = u238_single_resonance();
7069 let true_density = 0.0005_f64;
7070 let true_anorm = 0.85_f64;
7071 let true_ba = 0.03_f64;
7072 let true_bb = -0.01_f64;
7073 let true_bc = 0.0_f64;
7074 let c = 5.98_f64;
7075 let lam_ob = 2000.0_f64;
7076 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7077 let (t_inner, _) = synthetic_transmission(&data, true_density, &energies);
7078 let t_out: Vec<f64> = t_inner
7079 .iter()
7080 .zip(energies.iter())
7081 .map(|(&ti, &e)| true_anorm * ti + true_ba + true_bb / e.sqrt() + true_bc * e.sqrt())
7082 .collect();
7083 let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
7084 let sample_counts: Vec<f64> = t_out.iter().map(|&ti| c * lam_ob * ti).collect();
7085
7086 let bg = BackgroundConfig {
7087 anorm_init: 1.0,
7088 back_a_init: 0.0,
7089 back_b_init: 0.0,
7090 back_c_init: 0.0,
7091 back_d_init: 0.01,
7092 back_f_init: 1.0,
7093 fit_anorm: true,
7094 fit_back_a: true,
7095 fit_back_b: true,
7096 fit_back_c: true,
7097 fit_back_d: false,
7098 fit_back_f: false,
7099 };
7100
7101 let config = UnifiedFitConfig::new(
7102 energies,
7103 vec![data],
7104 vec!["U-238".into()],
7105 0.0,
7106 None,
7107 vec![0.001],
7108 )
7109 .unwrap()
7110 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7111 .with_counts_background(CountsBackgroundConfig {
7112 c,
7113 ..Default::default()
7114 })
7115 .with_transmission_background(bg)
7116 .with_counts_enable_polish(Some(true));
7124
7125 let input = InputData::Counts {
7126 sample_counts,
7127 open_beam_counts,
7128 };
7129 let r = fit_spectrum_typed(&input, &config).unwrap();
7130
7131 let dpd = r.deviance_per_dof.expect("joint-Poisson must report D/dof");
7141 assert!(
7142 dpd < 1.0,
7143 "D/dof = {dpd} unexpectedly large on noise-free fit — bg params not reaching objective?"
7144 );
7145 assert!(r.densities[0] > 1e-5, "density railed: {}", r.densities[0]);
7147 assert!(
7149 (r.anorm - 1.0).abs() > 0.05,
7150 "A_n did not move from init 1.0 (fitted={})",
7151 r.anorm
7152 );
7153 let bg_moved = r.background.iter().any(|v| v.abs() > 1e-4);
7155 assert!(
7156 bg_moved,
7157 "no bg parameter moved from init 0: {:?}",
7158 r.background
7159 );
7160 }
7161
7162 #[test]
7164 fn test_joint_poisson_requires_back_a_when_back_b_enabled() {
7165 let data = u238_single_resonance();
7166 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7167 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
7168 let ob: Vec<f64> = vec![500.0; energies.len()];
7169 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
7170
7171 let bg = BackgroundConfig {
7172 fit_anorm: true,
7174 fit_back_a: false,
7175 fit_back_b: true,
7176 fit_back_c: false,
7177 fit_back_d: false,
7178 fit_back_f: false,
7179 ..BackgroundConfig::default()
7180 };
7181
7182 let config = UnifiedFitConfig::new(
7183 energies,
7184 vec![data],
7185 vec!["U-238".into()],
7186 0.0,
7187 None,
7188 vec![0.001],
7189 )
7190 .unwrap()
7191 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7192 .with_counts_background(CountsBackgroundConfig {
7193 c: 1.0,
7194 ..Default::default()
7195 })
7196 .with_transmission_background(bg);
7197
7198 let input = InputData::Counts {
7199 sample_counts: s,
7200 open_beam_counts: ob,
7201 };
7202 let err = fit_spectrum_typed(&input, &config).unwrap_err();
7203 let msg = err.to_string();
7204 assert!(
7205 msg.contains("B_A"),
7206 "expected B_A pairing-rule rejection message, got: {msg}"
7207 );
7208 }
7209
7210 #[test]
7218 fn test_joint_poisson_uses_the_detector_background() {
7219 let data = u238_single_resonance();
7220 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7221 let true_density = 0.001;
7222 let (t, _) = synthetic_transmission(&data, true_density, &energies);
7223 let b_det = 60.0_f64;
7224 let flux: Vec<f64> = vec![2000.0 + b_det; energies.len()];
7226 let s: Vec<f64> = t.iter().map(|&ti| 2000.0 * ti + b_det).collect();
7227
7228 let fit = |background: Vec<f64>| {
7229 let config = UnifiedFitConfig::new(
7230 energies.clone(),
7231 vec![data.clone()],
7232 vec!["U-238".into()],
7233 0.0,
7234 None,
7235 vec![0.0005],
7236 )
7237 .unwrap()
7238 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7239 .with_counts_background(CountsBackgroundConfig {
7240 c: 1.0,
7241 ..Default::default()
7242 });
7243 let input = InputData::CountsWithNuisance {
7244 sample_counts: s.clone(),
7245 flux: flux.clone(),
7246 background,
7247 };
7248 let result = fit_spectrum_typed(&input, &config).expect("counts-KL fit runs");
7249 (result.densities[0] - true_density).abs() / true_density
7250 };
7251
7252 let declared = fit(vec![b_det; energies.len()]);
7253 let ignored = fit(vec![0.0; energies.len()]);
7254
7255 assert!(
7256 declared < 0.02,
7257 "declaring the background must recover density; bias {:.3} %",
7258 100.0 * declared
7259 );
7260 assert!(
7261 ignored > 4.0 * declared,
7262 "ignoring a background that IS in the counts gave bias {:.3} % against \
7263 {:.3} % when declared — the background is not reaching the likelihood",
7264 100.0 * ignored,
7265 100.0 * declared
7266 );
7267 }
7268
7269 #[test]
7271 fn test_joint_poisson_rejects_malformed_detector_background() {
7272 let data = u238_single_resonance();
7273 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7274 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
7275 let flux: Vec<f64> = vec![500.0; energies.len()];
7276 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
7277
7278 for bad in [f64::NAN, f64::INFINITY, -1.0] {
7279 let config = UnifiedFitConfig::new(
7280 energies.clone(),
7281 vec![data.clone()],
7282 vec!["U-238".into()],
7283 0.0,
7284 None,
7285 vec![0.001],
7286 )
7287 .unwrap()
7288 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7289 .with_counts_background(CountsBackgroundConfig {
7290 c: 1.0,
7291 ..Default::default()
7292 });
7293 let input = InputData::CountsWithNuisance {
7294 sample_counts: s.clone(),
7295 flux: flux.clone(),
7296 background: vec![bad; energies.len()],
7297 };
7298 let msg = fit_spectrum_typed(&input, &config).unwrap_err().to_string();
7299 assert!(
7300 msg.contains("finite and non-negative"),
7301 "expected a malformed-background rejection for {bad}, got: {msg}"
7302 );
7303 }
7304 }
7305
7306 #[test]
7308 fn test_joint_poisson_rejects_back_d_f() {
7309 let data = u238_single_resonance();
7310 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7311 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
7312 let ob: Vec<f64> = vec![500.0; energies.len()];
7313 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
7314
7315 let bg = BackgroundConfig {
7316 fit_back_d: true,
7317 ..BackgroundConfig::default()
7318 };
7319
7320 let config = UnifiedFitConfig::new(
7321 energies,
7322 vec![data],
7323 vec!["U-238".into()],
7324 0.0,
7325 None,
7326 vec![0.001],
7327 )
7328 .unwrap()
7329 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7330 .with_counts_background(CountsBackgroundConfig {
7331 c: 1.0,
7332 ..Default::default()
7333 })
7334 .with_transmission_background(bg);
7335
7336 let input = InputData::Counts {
7337 sample_counts: s,
7338 open_beam_counts: ob,
7339 };
7340 let err = fit_spectrum_typed(&input, &config).unwrap_err();
7341 assert!(
7342 err.to_string().contains("BackD"),
7343 "expected BackD/BackF rejection, got: {err}"
7344 );
7345 }
7346
7347 #[test]
7355 fn test_lm_transmission_rejects_partial_back_d_only() {
7356 let data = u238_single_resonance();
7357 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7358 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
7359
7360 let bg = BackgroundConfig {
7361 fit_back_d: true,
7364 fit_back_f: false,
7365 ..BackgroundConfig::default()
7366 };
7367
7368 let config = UnifiedFitConfig::new(
7369 energies,
7370 vec![data],
7371 vec!["U-238".into()],
7372 0.0,
7373 None,
7374 vec![0.001],
7375 )
7376 .unwrap()
7377 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7378 .with_transmission_background(bg);
7379
7380 let input = InputData::Transmission {
7381 transmission: t,
7382 uncertainty: u,
7383 };
7384 let err = fit_spectrum_typed(&input, &config).unwrap_err();
7385 let msg = err.to_string();
7386 assert!(
7387 msg.contains("fit_back_d") && msg.contains("fit_back_f"),
7388 "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
7389 );
7390 }
7391
7392 #[test]
7394 fn test_lm_transmission_rejects_partial_back_f_only() {
7395 let data = u238_single_resonance();
7396 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
7397 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
7398
7399 let bg = BackgroundConfig {
7400 fit_back_d: false,
7401 fit_back_f: true,
7402 ..BackgroundConfig::default()
7403 };
7404
7405 let config = UnifiedFitConfig::new(
7406 energies,
7407 vec![data],
7408 vec!["U-238".into()],
7409 0.0,
7410 None,
7411 vec![0.001],
7412 )
7413 .unwrap()
7414 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7415 .with_transmission_background(bg);
7416
7417 let input = InputData::Transmission {
7418 transmission: t,
7419 uncertainty: u,
7420 };
7421 let err = fit_spectrum_typed(&input, &config).unwrap_err();
7422 let msg = err.to_string();
7423 assert!(
7424 msg.contains("fit_back_d") && msg.contains("fit_back_f"),
7425 "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
7426 );
7427 }
7428
7429 #[test]
7435 fn test_unified_fit_config_fit_energy_range_round_trips() {
7436 let data = u238_single_resonance();
7437 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
7438 let cfg = UnifiedFitConfig::new(
7439 energies,
7440 vec![data],
7441 vec!["U-238".into()],
7442 0.0,
7443 None,
7444 vec![0.001],
7445 )
7446 .unwrap();
7447 assert_eq!(cfg.fit_energy_range(), None);
7448
7449 let cfg = cfg.with_fit_energy_range(Some((5.0, 50.0))).unwrap();
7450 assert_eq!(cfg.fit_energy_range(), Some((5.0, 50.0)));
7451
7452 let cfg = cfg.with_fit_energy_range(None).unwrap();
7454 assert_eq!(cfg.fit_energy_range(), None);
7455 }
7456
7457 #[test]
7461 fn test_unified_fit_config_fit_energy_range_rejects_invalid() {
7462 let data = u238_single_resonance();
7463 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
7464 let cfg = UnifiedFitConfig::new(
7465 energies,
7466 vec![data],
7467 vec!["U-238".into()],
7468 0.0,
7469 None,
7470 vec![0.001],
7471 )
7472 .unwrap();
7473
7474 let err = cfg
7476 .clone()
7477 .with_fit_energy_range(Some((10.0, 5.0)))
7478 .unwrap_err();
7479 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
7480
7481 let err = cfg
7483 .clone()
7484 .with_fit_energy_range(Some((5.0, 5.0)))
7485 .unwrap_err();
7486 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
7487
7488 let err = cfg
7490 .clone()
7491 .with_fit_energy_range(Some((f64::NAN, 5.0)))
7492 .unwrap_err();
7493 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
7494
7495 let err = cfg
7496 .clone()
7497 .with_fit_energy_range(Some((5.0, f64::INFINITY)))
7498 .unwrap_err();
7499 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
7500 }
7501
7502 #[test]
7507 fn test_fit_energy_range_lm_matches_subset_when_outside_negligible() {
7508 let data = u238_single_resonance();
7509 let true_density = 0.002;
7510 let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
7515 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
7516
7517 let cfg_full = UnifiedFitConfig::new(
7518 energies.clone(),
7519 vec![data.clone()],
7520 vec!["U-238".into()],
7521 0.0,
7522 None,
7523 vec![0.001],
7524 )
7525 .unwrap()
7526 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
7527
7528 let cfg_masked = UnifiedFitConfig::new(
7529 energies.clone(),
7530 vec![data],
7531 vec!["U-238".into()],
7532 0.0,
7533 None,
7534 vec![0.001],
7535 )
7536 .unwrap()
7537 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7538 .with_fit_energy_range(Some((4.0, 8.0)))
7539 .unwrap();
7540
7541 let input = InputData::Transmission {
7542 transmission: t,
7543 uncertainty: sigma,
7544 };
7545 let r_full = fit_spectrum_typed(&input, &cfg_full).unwrap();
7546 let r_masked = fit_spectrum_typed(&input, &cfg_masked).unwrap();
7547 assert!(r_full.converged && r_masked.converged);
7548
7549 let d_full = r_full.densities[0];
7550 let d_masked = r_masked.densities[0];
7551 let rel_err = (d_full - d_masked).abs() / d_full.abs();
7552 assert!(
7553 rel_err < 0.01,
7554 "fit_energy_range LM density {d_masked} should match full-grid {d_full} \
7555 within 1% (got rel_err = {rel_err})"
7556 );
7557 }
7558
7559 #[test]
7566 fn test_fit_energy_range_does_not_enable_transmission_kl() {
7567 let data = u238_single_resonance();
7568 let true_density = 0.002;
7569 let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
7570 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
7571
7572 let cfg = UnifiedFitConfig::new(
7573 energies,
7574 vec![data],
7575 vec!["U-238".into()],
7576 0.0,
7577 None,
7578 vec![0.001],
7579 )
7580 .unwrap()
7581 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7582 .with_fit_energy_range(Some((4.0, 8.0)))
7583 .unwrap();
7584
7585 let input = InputData::Transmission {
7586 transmission: t,
7587 uncertainty: sigma,
7588 };
7589 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
7590 let msg = err.to_string();
7591 assert!(
7592 msg.contains("normalized transmission") && msg.contains("Poisson"),
7593 "transmission plus Poisson/KL must be rejected as a route regardless \
7594 of fit_energy_range; got: {msg}"
7595 );
7596 }
7597
7598 #[test]
7603 fn test_fit_energy_range_lm_rejects_too_narrow() {
7604 let data = u238_single_resonance();
7605 let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
7609 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
7610 let cfg = UnifiedFitConfig::new(
7611 energies,
7612 vec![data],
7613 vec!["U-238".into()],
7614 0.0,
7615 None,
7616 vec![0.001],
7617 )
7618 .unwrap()
7619 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7620 .with_fit_energy_range(Some((4.6, 4.7)))
7621 .unwrap();
7622 let input = InputData::Transmission {
7623 transmission: t,
7624 uncertainty: sigma,
7625 };
7626 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
7627 let msg = err.to_string();
7628 assert!(
7629 msg.contains("fit_energy_range") && msg.contains("active bin"),
7630 "expected too-narrow-range rejection; got: {msg}"
7631 );
7632 }
7633
7634 #[test]
7637 fn test_fit_energy_range_jp_rejects_too_narrow() {
7638 let data = u238_single_resonance();
7639 let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
7640 let (sample, open_beam) = synthetic_counts(&data, 0.002, &energies, 1000.0);
7641 let cfg = UnifiedFitConfig::new(
7642 energies,
7643 vec![data],
7644 vec!["U-238".into()],
7645 0.0,
7646 None,
7647 vec![0.001],
7648 )
7649 .unwrap()
7650 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7651 .with_fit_energy_range(Some((4.6, 4.7)))
7652 .unwrap();
7653 let input = InputData::Counts {
7654 sample_counts: sample,
7655 open_beam_counts: open_beam,
7656 };
7657 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
7658 let msg = err.to_string();
7659 assert!(
7660 msg.contains("fit_energy_range") && msg.contains("active bin"),
7661 "expected too-narrow-range rejection; got: {msg}"
7662 );
7663 }
7664
7665 #[test]
7673 fn evaluate_jacobian_and_fisher_rejects_resolution() {
7674 use nereids_physics::resolution::{ResolutionFunction, ResolutionParams};
7675 let data = u238_single_resonance();
7676 let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
7677 let n_e = energies.len();
7678 let config = UnifiedFitConfig::new(
7679 energies,
7680 vec![data],
7681 vec!["U-238".into()],
7682 300.0,
7683 Some(ResolutionFunction::Gaussian(
7684 ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap(),
7685 )),
7686 vec![0.001],
7687 )
7688 .unwrap();
7689 let flux = vec![5000.0; n_e];
7690 let background = vec![10.0; n_e];
7691 let err = match evaluate_jacobian_and_fisher(&config, &flux, &background) {
7692 Ok(_) => panic!("counts-space Fisher helper accepted instrument resolution"),
7693 Err(err) => err,
7694 };
7695 let msg = err.to_string();
7696 assert!(
7697 msg.contains("instrument resolution") && msg.contains("separate-arm model"),
7698 "expected physical counts-response rejection, got: {msg}"
7699 );
7700 }
7701
7702 #[test]
7712 fn fit_spectrum_typed_energy_scale_lm_recovers_calibration() {
7713 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
7715 let true_density = 0.05_f64;
7716 let (t0_true, ls_true) = (1.5_f64, 1.004_f64);
7717 let model = EnergyScaleTransmissionModel::new(
7720 Arc::new(vec![data.clone()]),
7721 Arc::new(vec![0]),
7722 Arc::new(vec![1.0]),
7723 293.6,
7724 energies.clone(),
7725 25.0,
7726 1,
7727 2,
7728 None,
7729 );
7730 let measured = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
7731 let uncertainty = vec![1e-3; energies.len()];
7732 let config = UnifiedFitConfig::new(
7733 energies,
7734 vec![data],
7735 vec!["Hf-178".into()],
7736 293.6,
7737 None,
7738 vec![0.04], )
7740 .unwrap()
7741 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7742 .with_energy_scale(0.0, 1.0, 25.0); let input = InputData::Transmission {
7744 transmission: measured,
7745 uncertainty,
7746 };
7747 let result = fit_spectrum_typed(&input, &config).unwrap();
7748 let t0 = result
7749 .t0_us
7750 .expect("t0_us populated when fit_energy_scale=true");
7751 let ls = result
7752 .l_scale
7753 .expect("l_scale populated when fit_energy_scale=true");
7754 assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
7755 assert!(result.converged, "energy-scale LM fit should converge");
7756 assert!(
7757 (result.densities[0] - true_density).abs() / true_density < 0.10,
7758 "density: fitted={}, true={true_density}",
7759 result.densities[0]
7760 );
7761 }
7762
7763 #[test]
7768 fn fit_spectrum_typed_energy_scale_counts_kl_seeds_via_proxy() {
7769 let data = hf178_mlbw_two_resonances();
7770 let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
7771 let true_density = 0.05_f64;
7772 let (t0_true, ls_true) = (1.0_f64, 1.003_f64);
7773 let model = EnergyScaleTransmissionModel::new(
7774 Arc::new(vec![data.clone()]),
7775 Arc::new(vec![0]),
7776 Arc::new(vec![1.0]),
7777 293.6,
7778 energies.clone(),
7779 25.0,
7780 1,
7781 2,
7782 None,
7783 );
7784 let t = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
7785 let flux: Vec<f64> = vec![5000.0; energies.len()];
7789 let background: Vec<f64> = vec![0.0; energies.len()];
7790 let sample: Vec<f64> = t
7791 .iter()
7792 .zip(flux.iter())
7793 .map(|(&ti, &fi)| fi * ti)
7794 .collect();
7795 let config = UnifiedFitConfig::new(
7796 energies,
7797 vec![data],
7798 vec!["Hf-178".into()],
7799 293.6,
7800 None,
7801 vec![0.04],
7802 )
7803 .unwrap()
7804 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7805 .with_energy_scale(0.0, 1.0, 25.0);
7806 let input = InputData::CountsWithNuisance {
7807 sample_counts: sample,
7808 flux,
7809 background,
7810 };
7811 let result = fit_spectrum_typed(&input, &config).unwrap();
7812 let t0 = result
7813 .t0_us
7814 .expect("t0_us populated when fit_energy_scale=true");
7815 let ls = result
7816 .l_scale
7817 .expect("l_scale populated when fit_energy_scale=true");
7818 assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
7819 }
7820
7821 #[test]
7822 fn validate_precomputed_cross_sections_error_branches() {
7823 use nereids_physics::resolution::{ResolutionFunction, ResolutionParams};
7824 let data = u238_single_resonance();
7825 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
7826 let n_e = energies.len();
7827 let base = UnifiedFitConfig::new(
7828 energies.clone(),
7829 vec![data.clone()],
7830 vec!["U-238".into()],
7831 0.0,
7832 None,
7833 vec![0.001],
7834 )
7835 .unwrap();
7836 let identity = || Arc::new(WorkingGridLayout::identity(&energies));
7837 let cfg =
7838 |base: &UnifiedFitConfig, sigma: Vec<Vec<f64>>, layout: Arc<WorkingGridLayout>| {
7839 base.clone().with_precomputed_cross_sections(PrecomputedXs {
7840 sigma: Arc::new(sigma),
7841 layout,
7842 })
7843 };
7844 let expect =
7845 |c: &UnifiedFitConfig, needle: &str| match validate_precomputed_cross_sections(c) {
7846 Err(PipelineError::ShapeMismatch(m)) => {
7847 assert!(m.contains(needle), "expected {needle:?}, got: {m}")
7848 }
7849 other => panic!("expected ShapeMismatch({needle:?}), got {other:?}"),
7850 };
7851 assert!(
7852 validate_precomputed_cross_sections(&cfg(&base, vec![vec![1.0; n_e]], identity()))
7853 .is_ok()
7854 );
7855 expect(&cfg(&base, vec![], identity()), "must not be empty");
7856 expect(
7857 &cfg(&base, vec![vec![1.0; n_e - 1]], identity()),
7858 "its grid has",
7859 );
7860 let mut nan_row = vec![1.0; n_e];
7861 nan_row[3] = f64::NAN;
7862 expect(&cfg(&base, vec![nan_row], identity()), "non-finite");
7863 expect(
7864 &cfg(&base, vec![vec![1.0; n_e], vec![1.0; n_e]], identity()),
7865 "rows but expected",
7866 );
7867 let extended = Arc::new(WorkingGridLayout {
7868 energies: {
7869 let mut e = energies.clone();
7870 e.insert(0, 0.95);
7871 e.push(2.05);
7872 e
7873 },
7874 data_indices: (1..=n_e).collect(),
7875 });
7876 expect(
7877 &cfg(&base, vec![vec![1.0; n_e + 2]], extended),
7878 "is not the working grid",
7879 );
7880
7881 let res =
7882 ResolutionFunction::Gaussian(ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap());
7883 let with_res = UnifiedFitConfig::new(
7884 energies.clone(),
7885 vec![data.clone()],
7886 vec!["U-238".into()],
7887 0.0,
7888 Some(res.clone()),
7889 vec![0.001],
7890 )
7891 .unwrap();
7892 expect(
7893 &cfg(&with_res, vec![vec![1.0; n_e]], identity()),
7894 "is not the working grid",
7895 );
7896 let working = Arc::new(
7897 phys_transmission::resolution_working_grid(
7898 &energies,
7899 Some(&InstrumentParams { resolution: res }),
7900 &[&data],
7901 )
7902 .unwrap(),
7903 );
7904 assert!(!working.is_identity());
7905 let n_work = working.energies.len();
7906 assert!(
7907 validate_precomputed_cross_sections(&cfg(
7908 &with_res,
7909 vec![vec![1.0; n_work]],
7910 Arc::clone(&working)
7911 ))
7912 .is_ok()
7913 );
7914 let with_indices = |data_indices: Vec<usize>| {
7915 Arc::new(WorkingGridLayout {
7916 energies: working.energies.clone(),
7917 data_indices,
7918 })
7919 };
7920 let mut shortened = working.data_indices.clone();
7921 shortened.pop();
7922 let mut out_of_range = working.data_indices.clone();
7923 out_of_range[0] = n_work + 5;
7924 let mut permuted = working.data_indices.clone();
7925 permuted.reverse();
7926 for bad in [shortened, out_of_range, permuted] {
7927 expect(
7928 &cfg(&with_res, vec![vec![1.0; n_work]], with_indices(bad)),
7929 "is not the working grid",
7930 );
7931 }
7932 }
7933
7934 #[test]
7935 fn evaluate_jacobian_and_fisher_identity_and_temperature_paths() {
7936 let data = u238_single_resonance();
7937 let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
7938 let n_e = energies.len();
7939 let flux = vec![5000.0; n_e];
7940 let background = vec![10.0; n_e];
7941 let cfg_id = UnifiedFitConfig::new(
7943 energies.clone(),
7944 vec![data.clone()],
7945 vec!["U-238".into()],
7946 300.0,
7947 None,
7948 vec![0.001],
7949 )
7950 .unwrap();
7951 let r_id = evaluate_jacobian_and_fisher(&cfg_id, &flux, &background).unwrap();
7952 assert_eq!(r_id.model_prediction.len(), n_e);
7953 let f00 = r_id.fisher.get(0, 0);
7954 assert!(
7955 f00.is_finite() && f00 > 0.0,
7956 "identity-path Fisher[0,0]={f00}"
7957 );
7958 let cfg_t = UnifiedFitConfig::new(
7959 energies,
7960 vec![data],
7961 vec!["U-238".into()],
7962 300.0,
7963 None,
7964 vec![0.001],
7965 )
7966 .unwrap()
7967 .with_fit_temperature(true);
7968 let r_t = evaluate_jacobian_and_fisher(&cfg_t, &flux, &background).unwrap();
7969 assert_eq!(
7970 r_t.param_names.len(),
7971 2,
7972 "density + temperature free params"
7973 );
7974 assert!(r_t.model_prediction.iter().all(|v| v.is_finite()));
7975 }
7976
7977 #[test]
7982 fn fit_spectrum_typed_grouped_precomputed_collapses_by_groups() {
7983 use nereids_core::types::{Isotope, IsotopeGroup};
7984 let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
7985 let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);
7986 let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
7987 let true_density = 0.0005_f64;
7988 let sample = phys_transmission::SampleParams::new(
7989 0.0,
7990 vec![
7991 (rd1.clone(), true_density * 0.6),
7992 (rd2.clone(), true_density * 0.4),
7993 ],
7994 )
7995 .unwrap();
7996 let transmission = phys_transmission::forward_model(&energies, &sample, None).unwrap();
7997 let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();
7998 let per_member = phys_transmission::broadened_cross_sections(
8000 &energies,
8001 &[rd1.clone(), rd2.clone()],
8002 0.0,
8003 None,
8004 None,
8005 )
8006 .unwrap();
8007 let iso1 = Isotope::new(92, 235).unwrap();
8008 let iso2 = Isotope::new(92, 238).unwrap();
8009 let group = IsotopeGroup::custom("U".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();
8010 let config = UnifiedFitConfig::new(
8011 energies,
8012 vec![rd1.clone()],
8013 vec!["placeholder".into()],
8014 0.0,
8015 None,
8016 vec![0.001],
8017 )
8018 .unwrap()
8019 .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
8020 .unwrap();
8021 let per_member = table_on_data_grid(&config, per_member);
8022 let config = config
8023 .with_precomputed_cross_sections(per_member)
8024 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
8025 let input = InputData::Transmission {
8026 transmission,
8027 uncertainty,
8028 };
8029 let result = fit_spectrum_typed(&input, &config).unwrap();
8030 assert_eq!(result.densities.len(), 1, "one group density");
8031 assert!(
8032 (result.densities[0] - true_density).abs() / true_density < 0.01,
8033 "grouped+precomputed density: fitted={}, true={true_density}",
8034 result.densities[0]
8035 );
8036 }
8037
8038 #[test]
8042 fn peak_match_energy_scale_seed_rejects_out_of_bounds() {
8043 let data = hf178_mlbw_two_resonances();
8044 let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
8045 let model = EnergyScaleTransmissionModel::new(
8046 Arc::new(vec![data.clone()]),
8047 Arc::new(vec![0]),
8048 Arc::new(vec![1.0]),
8049 293.6,
8050 energies.clone(),
8051 25.0,
8052 1,
8053 2,
8054 None,
8055 );
8056 let t_obs = model.evaluate(&[0.05, 0.0, 1.03]).unwrap();
8058 let config = UnifiedFitConfig::new(
8059 energies.clone(),
8060 vec![data],
8061 vec!["Hf-178".into()],
8062 293.6,
8063 None,
8064 vec![0.05],
8065 )
8066 .unwrap()
8067 .with_energy_scale(0.0, 1.0, 25.0);
8068 assert!(
8069 peak_match_energy_scale_seed(
8070 &t_obs,
8071 config.energies(),
8072 &config,
8073 25.0,
8074 (-10.0, 10.0),
8075 (0.99, 1.01),
8076 )
8077 .is_none(),
8078 "an out-of-bounds calibration fit must be rejected (cold-start fallback)"
8079 );
8080 }
8081
8082 #[test]
8088 fn test_count_free_params_drops_fixed_densities() {
8089 let data = u238_single_resonance();
8090 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
8091 let base = UnifiedFitConfig::new(
8092 energies,
8093 vec![data.clone(), data],
8094 vec!["a".into(), "b".into()],
8095 300.0,
8096 None,
8097 vec![0.001, 0.001],
8098 )
8099 .unwrap()
8100 .with_fit_temperature(true);
8101
8102 assert_eq!(count_free_params(&base), 3);
8104 assert_eq!(count_free_params(&base.clone().with_fix_densities(true)), 1);
8106 let one_fixed = base.with_density_free(vec![false, true]).unwrap();
8108 assert_eq!(count_free_params(&one_fixed), 2);
8109 }
8110
8111 #[test]
8114 fn test_with_density_free_validates_and_normalises() {
8115 let data = u238_single_resonance();
8116 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
8117 let cfg = UnifiedFitConfig::new(
8118 energies,
8119 vec![data],
8120 vec!["a".into()],
8121 300.0,
8122 None,
8123 vec![0.001],
8124 )
8125 .unwrap();
8126 assert!(matches!(
8128 cfg.clone().with_density_free(vec![true, false]),
8129 Err(FitConfigError::DensityCountMismatch { .. })
8130 ));
8131 let all_free = cfg.with_density_free(vec![true]).unwrap();
8133 assert_eq!(all_free.n_free_density_params(), 1);
8134 assert!(!all_free.density_is_fixed(0));
8135 }
8136
8137 #[test]
8143 fn test_with_groups_rejects_prior_density_freeze() {
8144 use nereids_core::types::{Isotope, IsotopeGroup};
8145 let data = u238_single_resonance();
8146 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
8147 let iso = Isotope::new(92, 238).unwrap();
8148 let group = IsotopeGroup::custom("U-238".into(), vec![(iso, 1.0)]).unwrap();
8149
8150 let cfg = UnifiedFitConfig::new(
8151 energies,
8152 vec![data.clone()],
8153 vec!["placeholder".into()],
8154 300.0,
8155 None,
8156 vec![0.001],
8157 )
8158 .unwrap()
8159 .with_fix_densities(true);
8161 assert!(cfg.density_is_fixed(0), "mask set before grouping");
8162
8163 assert!(matches!(
8165 cfg.with_groups(&[(&group, &[data])], vec![0.001]),
8166 Err(FitConfigError::DensityFreezeBeforeGroups)
8167 ));
8168 }
8169
8170 #[test]
8173 fn test_freeze_after_groups_succeeds() {
8174 use nereids_core::types::{Isotope, IsotopeGroup};
8175 let data = u238_single_resonance();
8176 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
8177 let iso = Isotope::new(92, 238).unwrap();
8178 let group = IsotopeGroup::custom("U-238".into(), vec![(iso, 1.0)]).unwrap();
8179
8180 let cfg = UnifiedFitConfig::new(
8181 energies,
8182 vec![data.clone()],
8183 vec!["placeholder".into()],
8184 300.0,
8185 None,
8186 vec![0.001],
8187 )
8188 .unwrap()
8189 .with_groups(&[(&group, &[data])], vec![0.001])
8190 .unwrap()
8191 .with_fix_densities(true);
8192 assert_eq!(cfg.n_density_params(), 1);
8193 assert!(
8194 cfg.density_is_fixed(0),
8195 "group density frozen when freeze follows grouping"
8196 );
8197 assert_eq!(cfg.n_free_density_params(), 0);
8198 }
8199
8200 #[test]
8204 fn test_all_frozen_no_free_param_rejected() {
8205 let data = u238_single_resonance();
8206 let true_density = 0.0005;
8207 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8208 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, 300.0, &energies);
8209 let config = UnifiedFitConfig::new(
8210 energies,
8211 vec![data],
8212 vec!["U-238".into()],
8213 300.0,
8214 None,
8215 vec![true_density],
8216 )
8217 .unwrap()
8218 .with_fix_densities(true);
8220 assert_eq!(count_free_params(&config), 0);
8221 let input = InputData::Transmission {
8222 transmission: t,
8223 uncertainty: sigma,
8224 };
8225 assert!(matches!(
8226 fit_spectrum_typed(&input, &config),
8227 Err(PipelineError::InvalidParameter(_))
8228 ));
8229 }
8230
8231 #[test]
8236 fn test_fix_densities_recovers_temperature_holds_density() {
8237 let data = u238_single_resonance();
8238 let true_density = 0.0005;
8239 let true_temp = 350.0;
8240 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8241 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
8242
8243 let config = UnifiedFitConfig::new(
8244 energies,
8245 vec![data],
8246 vec!["U-238".into()],
8247 300.0, None,
8249 vec![true_density], )
8251 .unwrap()
8252 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
8253 .with_fit_temperature(true)
8254 .with_fix_densities(true);
8255
8256 let input = InputData::Transmission {
8257 transmission: t,
8258 uncertainty: sigma,
8259 };
8260 let result = fit_spectrum_typed(&input, &config).unwrap();
8261 assert!(result.converged, "T-only fit should converge");
8262
8263 assert_eq!(
8265 result.densities[0], true_density,
8266 "frozen density must not move"
8267 );
8268 let fitted_temp = result
8270 .temperature_k
8271 .expect("temperature_k should be Some when fit_temperature=true");
8272 assert!(
8273 (fitted_temp - true_temp).abs() < 1.0,
8274 "temperature: fitted={fitted_temp}, true={true_temp}"
8275 );
8276
8277 let t_unc = result
8283 .temperature_k_unc
8284 .expect("temperature_k_unc should be Some for a converged T fit");
8285 assert!(
8286 t_unc.is_finite() && t_unc > 0.0,
8287 "frozen-density thermometry must report a finite positive temperature σ, got {t_unc}"
8288 );
8289 let dens_unc = result
8290 .uncertainties
8291 .as_ref()
8292 .expect("converged fit has density uncertainties");
8293 assert!(
8294 dens_unc[0].is_nan(),
8295 "frozen density must report NaN σ (no covariance column), got {}",
8296 dens_unc[0]
8297 );
8298 }
8299
8300 #[test]
8303 fn test_fix_densities_does_not_enable_transmission_kl() {
8304 let data = u238_single_resonance();
8305 let true_density = 0.0005;
8306 let true_temp = 350.0;
8307 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8308 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
8309
8310 let config = UnifiedFitConfig::new(
8311 energies,
8312 vec![data],
8313 vec!["U-238".into()],
8314 300.0,
8315 None,
8316 vec![true_density],
8317 )
8318 .unwrap()
8319 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
8320 .with_fit_temperature(true)
8321 .with_fix_densities(true);
8322
8323 let input = InputData::Transmission {
8324 transmission: t,
8325 uncertainty: sigma,
8326 };
8327 let error = fit_spectrum_typed(&input, &config)
8328 .expect_err("transmission plus Poisson/KL must be rejected");
8329 assert!(error.to_string().contains("normalized transmission"));
8330 }
8331
8332 #[test]
8339 fn test_extract_result_maps_uncertainties_past_frozen_density() {
8340 let data = u238_single_resonance();
8341 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
8342 let config = UnifiedFitConfig::new(
8344 energies,
8345 vec![data.clone(), data],
8346 vec!["a".into(), "b".into()],
8347 300.0,
8348 None,
8349 vec![0.001, 0.002],
8350 )
8351 .unwrap()
8352 .with_fit_temperature(true)
8353 .with_density_free(vec![false, true])
8354 .unwrap();
8355
8356 let result = LmResult {
8359 chi_squared: 1.0,
8360 reduced_chi_squared: 1.0,
8361 iterations: 5,
8362 converged: true,
8363 params: vec![0.001, 0.002, 350.0],
8364 covariance: Some(lm::FlatMatrix::zeros(2, 2)),
8365 uncertainties: Some(vec![0.02, 4.0]),
8366 };
8367
8368 let extracted = extract_result(&config, &result, 2, &[1, 2], None, None).unwrap();
8369 let unc = extracted
8370 .uncertainties
8371 .expect("converged fit surfaces density uncertainties");
8372 assert!(
8373 unc[0].is_nan(),
8374 "frozen density d0 must report NaN σ, got {}",
8375 unc[0]
8376 );
8377 assert_eq!(unc[1], 0.02, "free density d1 σ must come from free slot 0");
8378 assert_eq!(
8379 extracted.temperature_k_unc,
8380 Some(4.0),
8381 "temperature σ must come from free slot 1, not the missing full index 2"
8382 );
8383 }
8384
8385 #[test]
8389 fn test_fix_densities_biased_density_gives_finite_documented_bias() {
8390 let data = u238_single_resonance();
8391 let true_density = 0.0005;
8392 let true_temp = 350.0;
8393 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8394 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
8395
8396 let config = UnifiedFitConfig::new(
8397 energies,
8398 vec![data],
8399 vec!["U-238".into()],
8400 300.0,
8401 None,
8402 vec![0.9 * true_density], )
8404 .unwrap()
8405 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
8406 .with_fit_temperature(true)
8407 .with_fix_densities(true);
8408
8409 let input = InputData::Transmission {
8410 transmission: t,
8411 uncertainty: sigma,
8412 };
8413 let result = fit_spectrum_typed(&input, &config).unwrap();
8414 let fitted_temp = result.temperature_k.expect("temperature_k should be Some");
8415 assert_eq!(result.densities[0], 0.9 * true_density);
8417 assert!(
8418 fitted_temp.is_finite() && fitted_temp > 1.0,
8419 "biased-density fit must yield a finite temperature, got {fitted_temp}"
8420 );
8421 }
8422
8423 const BL_TRUE: [f64; 3] = [1.02, -0.03, 0.01];
8428
8429 fn baseline_at(e: f64, e_ref: f64, b: &[f64; 3]) -> f64 {
8430 let z = (e / e_ref).ln();
8431 b[0] + b[1] * z + b[2] * z * z
8432 }
8433
8434 fn apply_truth_baseline(t: &[f64], energies: &[f64]) -> Vec<f64> {
8438 let e_ref = nereids_fitting::transmission_model::baseline_reference_energy(energies);
8439 let out: Vec<f64> = t
8440 .iter()
8441 .zip(energies.iter())
8442 .map(|(&ti, &e)| ti * baseline_at(e, e_ref, &BL_TRUE))
8443 .collect();
8444 let max_rel = t
8445 .iter()
8446 .zip(out.iter())
8447 .map(|(&a, &b)| ((b - a) / a.max(1e-12)).abs())
8448 .fold(0.0_f64, f64::max);
8449 assert!(
8450 max_rel > 0.01,
8451 "non-vacuity pre-check: injected baseline moved the data by only \
8452 {max_rel:.2e} (max relative) — the closed loop would not \
8453 distinguish baseline-on from baseline-off"
8454 );
8455 out
8456 }
8457
8458 #[test]
8459 fn baseline_lm_closed_loop_recovers_temperature_and_coefficients() {
8460 let data = u238_single_resonance();
8461 let true_density = 0.002;
8462 let true_temp = 600.0;
8463 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8464 let (t_pure, _) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
8465 let measured = apply_truth_baseline(&t_pure, &energies);
8466 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
8467 let e_ref = nereids_fitting::transmission_model::baseline_reference_energy(&energies);
8468
8469 let config = UnifiedFitConfig::new(
8473 energies,
8474 vec![data],
8475 vec!["U-238".into()],
8476 500.0,
8477 None,
8478 vec![true_density],
8479 )
8480 .unwrap()
8481 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
8482 .with_fit_temperature(true)
8483 .with_fix_densities(true)
8484 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8485
8486 let input = InputData::Transmission {
8487 transmission: measured,
8488 uncertainty: sigma,
8489 };
8490 let result = fit_spectrum_typed(&input, &config).unwrap();
8491 assert!(result.converged, "baseline LM closed loop should converge");
8492
8493 let fitted_temp = result.temperature_k.expect("temperature fitted");
8494 assert!(
8495 (fitted_temp - true_temp).abs() < 10.0,
8496 "temperature: fitted={fitted_temp}, true={true_temp}"
8497 );
8498
8499 let b = result.baseline.expect("baseline fitted");
8500 for (i, (&fitted, &truth)) in b.iter().zip(BL_TRUE.iter()).enumerate() {
8501 assert!(
8502 (fitted - truth).abs() < 1e-2,
8503 "b{i}: fitted={fitted}, true={truth}"
8504 );
8505 }
8506 let bl_cfg = MultiplicativeBaselineConfig::default();
8507 for (i, (&fitted, bounds)) in b
8508 .iter()
8509 .zip([bl_cfg.b0_bounds, bl_cfg.b1_bounds, bl_cfg.b2_bounds].iter())
8510 .enumerate()
8511 {
8512 assert!(
8513 fitted >= bounds.0 && fitted <= bounds.1,
8514 "b{i} = {fitted} escaped its bounds {bounds:?}"
8515 );
8516 }
8517 let reported_e_ref = result.baseline_e_ref_ev.expect("e_ref reported");
8518 assert!(
8519 (reported_e_ref - e_ref).abs() < 1e-12,
8520 "reported E_ref {reported_e_ref} != geometric midpoint {e_ref}"
8521 );
8522 assert!(
8523 result.warnings.is_empty(),
8524 "no degenerate trio here — warnings should be empty, got {:?}",
8525 result.warnings
8526 );
8527 }
8528
8529 #[test]
8530 fn baseline_does_not_enable_transmission_kl() {
8531 let data = u238_single_resonance();
8532 let true_density = 0.002;
8533 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8534 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
8535 let measured = apply_truth_baseline(&t_pure, &energies);
8536 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
8537
8538 let config = UnifiedFitConfig::new(
8539 energies,
8540 vec![data],
8541 vec!["U-238".into()],
8542 0.0,
8543 None,
8544 vec![0.001],
8545 )
8546 .unwrap()
8547 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
8548 max_iter: 500,
8549 ..PoissonConfig::default()
8550 }))
8551 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8552
8553 let input = InputData::Transmission {
8554 transmission: measured,
8555 uncertainty: sigma,
8556 };
8557 let error = fit_spectrum_typed(&input, &config)
8558 .expect_err("a baseline must not enable transmission plus Poisson/KL");
8559 assert!(error.to_string().contains("normalized transmission"));
8560 }
8561
8562 #[test]
8563 fn baseline_counts_jp_closed_loop() {
8564 let data = u238_single_resonance();
8565 let true_density = 0.002;
8566 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8567 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
8568 let t_baselined = apply_truth_baseline(&t_pure, &energies);
8569 let i0 = 10_000.0;
8570 let open_beam: Vec<f64> = vec![i0; t_baselined.len()];
8571 let sample: Vec<f64> = t_baselined
8572 .iter()
8573 .map(|&v| (v * i0).round().max(0.0))
8574 .collect();
8575
8576 let config = UnifiedFitConfig::new(
8577 energies,
8578 vec![data],
8579 vec!["U-238".into()],
8580 0.0,
8581 None,
8582 vec![0.001],
8583 )
8584 .unwrap()
8585 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
8586 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8587
8588 let input = InputData::Counts {
8589 sample_counts: sample,
8590 open_beam_counts: open_beam,
8591 };
8592 let result = fit_spectrum_typed(&input, &config).unwrap();
8593 assert!(result.converged, "baseline JP closed loop should converge");
8594 assert!(
8595 result.deviance_per_dof.is_some(),
8596 "counts-KL path reports deviance"
8597 );
8598 let fitted = result.densities[0];
8599 assert!(
8600 (fitted - true_density).abs() / true_density < 0.02,
8601 "density: fitted={fitted}, true={true_density}"
8602 );
8603 let b = result.baseline.expect("baseline fitted");
8604 for (i, (&fit_b, &truth)) in b.iter().zip(BL_TRUE.iter()).enumerate() {
8605 assert!(
8606 (fit_b - truth).abs() < 1e-2,
8607 "b{i}: fitted={fit_b}, true={truth}"
8608 );
8609 }
8610 }
8611
8612 #[test]
8616 fn baseline_rejects_free_anorm_on_all_paths() {
8617 let data = u238_single_resonance();
8618 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
8619 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
8620
8621 let base_config = UnifiedFitConfig::new(
8622 energies.clone(),
8623 vec![data],
8624 vec!["U-238".into()],
8625 0.0,
8626 None,
8627 vec![0.001],
8628 )
8629 .unwrap()
8630 .with_transmission_background(BackgroundConfig::default())
8632 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8633
8634 let transmission_input = InputData::Transmission {
8635 transmission: t.clone(),
8636 uncertainty: sigma,
8637 };
8638 let counts_input = InputData::Counts {
8639 sample_counts: t.iter().map(|&v| (v * 1000.0).round()).collect(),
8640 open_beam_counts: vec![1000.0; energies.len()],
8641 };
8642
8643 for (label, input, solver) in [
8644 (
8645 "LM transmission",
8646 &transmission_input,
8647 SolverConfig::LevenbergMarquardt(LmConfig::default()),
8648 ),
8649 (
8650 "joint-Poisson counts",
8651 &counts_input,
8652 SolverConfig::PoissonKL(PoissonConfig::default()),
8653 ),
8654 ] {
8655 let config = base_config.clone().with_solver(solver);
8656 let err = fit_spectrum_typed(input, &config)
8657 .expect_err(&format!("{label}: free Anorm + baseline must be rejected"));
8658 let msg = err.to_string();
8659 assert!(
8660 msg.contains("Anorm") && msg.contains("fit_anorm = false"),
8661 "{label}: rejection must name the degeneracy and the fix, got: {msg}"
8662 );
8663 }
8664 }
8665
8666 #[test]
8670 fn baseline_with_fixed_anorm_abc_accepted() {
8671 let data = u238_single_resonance();
8672 let true_density = 0.002;
8673 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8674 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
8675 let measured = apply_truth_baseline(&t_pure, &energies);
8676 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
8677
8678 let config = UnifiedFitConfig::new(
8679 energies,
8680 vec![data],
8681 vec!["U-238".into()],
8682 0.0,
8683 None,
8684 vec![0.001],
8685 )
8686 .unwrap()
8687 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
8688 .with_transmission_background(BackgroundConfig {
8689 fit_anorm: false,
8690 ..BackgroundConfig::default()
8691 })
8692 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8693
8694 let input = InputData::Transmission {
8695 transmission: measured,
8696 uncertainty: sigma,
8697 };
8698 let result = fit_spectrum_typed(&input, &config).unwrap();
8699 assert!(result.converged, "fixed-Anorm + ABC + baseline should fit");
8700 assert!(
8701 (result.anorm - 1.0).abs() < f64::EPSILON,
8702 "Anorm was frozen at 1.0, got {}",
8703 result.anorm
8704 );
8705 assert!(result.baseline.is_some(), "baseline block reported");
8706 let b = result.baseline.unwrap();
8709 assert!(
8710 (b[0] - BL_TRUE[0]).abs() < 0.02,
8711 "b0: fitted={}, true={}",
8712 b[0],
8713 BL_TRUE[0]
8714 );
8715 }
8716
8717 #[test]
8718 fn baseline_validation_rejects_bad_configs() {
8719 let data = u238_single_resonance();
8720 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
8721 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
8722 let input = InputData::Transmission {
8723 transmission: t,
8724 uncertainty: sigma,
8725 };
8726 let mk_config = |bl: MultiplicativeBaselineConfig| {
8727 UnifiedFitConfig::new(
8728 energies.clone(),
8729 vec![data.clone()],
8730 vec!["U-238".into()],
8731 0.0,
8732 None,
8733 vec![0.001],
8734 )
8735 .unwrap()
8736 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
8737 .with_multiplicative_baseline(bl)
8738 };
8739
8740 let err = fit_spectrum_typed(
8742 &input,
8743 &mk_config(MultiplicativeBaselineConfig {
8744 b0_init: f64::NAN,
8745 ..MultiplicativeBaselineConfig::default()
8746 }),
8747 )
8748 .expect_err("NaN b0_init must be rejected");
8749 assert!(err.to_string().contains("b0_init must be finite"));
8750
8751 let err = fit_spectrum_typed(
8753 &input,
8754 &mk_config(MultiplicativeBaselineConfig {
8755 b1_bounds: (0.05, -0.05),
8756 ..MultiplicativeBaselineConfig::default()
8757 }),
8758 )
8759 .expect_err("reversed b1_bounds must be rejected");
8760 assert!(err.to_string().contains("b1_bounds"));
8761
8762 let err = fit_spectrum_typed(
8764 &input,
8765 &mk_config(MultiplicativeBaselineConfig {
8766 b0_init: 1.5,
8767 ..MultiplicativeBaselineConfig::default()
8768 }),
8769 )
8770 .expect_err("out-of-bounds b0_init must be rejected");
8771 assert!(err.to_string().contains("outside"));
8772
8773 let wide_energies: Vec<f64> = (0..51).map(|i| 1.0 * 1.1f64.powi(i)).collect();
8778 let (wt, ws) = synthetic_transmission(&data, 0.002, &wide_energies);
8779 let wide_input = InputData::Transmission {
8780 transmission: wt,
8781 uncertainty: ws,
8782 };
8783 let wide_config = UnifiedFitConfig::new(
8784 wide_energies,
8785 vec![data.clone()],
8786 vec!["U-238".into()],
8787 0.0,
8788 None,
8789 vec![0.001],
8790 )
8791 .unwrap()
8792 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
8793 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
8794 b0_init: 0.3,
8795 b1_init: 0.4,
8796 b0_bounds: (0.1, 1.1),
8797 b1_bounds: (-1.0, 1.0),
8798 ..MultiplicativeBaselineConfig::default()
8799 });
8800 let err = fit_spectrum_typed(&wide_input, &wide_config)
8801 .expect_err("non-positive initial B(E) must be rejected");
8802 assert!(err.to_string().contains("not strictly"), "got: {err}");
8803 }
8804
8805 #[test]
8806 fn degenerate_trio_produces_warning() {
8807 let data = u238_single_resonance();
8808 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
8809
8810 let trio = UnifiedFitConfig::new(
8811 energies.clone(),
8812 vec![data.clone()],
8813 vec!["U-238".into()],
8814 300.0,
8815 None,
8816 vec![0.001],
8817 )
8818 .unwrap()
8819 .with_fit_temperature(true)
8820 .with_transmission_background(BackgroundConfig::default());
8821 let w = degenerate_normalization_warning(&trio)
8822 .expect("free Anorm + free T + free density must warn");
8823 assert!(w.contains("degenerate"), "warning names the failure: {w}");
8824
8825 let frozen_density = trio.clone().with_fix_densities(true);
8827 assert!(degenerate_normalization_warning(&frozen_density).is_none());
8828 let no_temp = trio.clone().with_fit_temperature(false);
8829 assert!(degenerate_normalization_warning(&no_temp).is_none());
8830 let fixed_anorm = trio.clone().with_transmission_background(BackgroundConfig {
8831 fit_anorm: false,
8832 ..BackgroundConfig::default()
8833 });
8834 assert!(degenerate_normalization_warning(&fixed_anorm).is_none());
8835
8836 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
8839 let input = InputData::Transmission {
8840 transmission: t,
8841 uncertainty: sigma,
8842 };
8843 let config = trio.with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
8844 max_iter: 5,
8845 ..LmConfig::default()
8846 }));
8847 let result = fit_spectrum_typed(&input, &config).unwrap();
8848 assert!(
8849 result.warnings.iter().any(|w| w.contains("degenerate")),
8850 "fit result must carry the degenerate-trio warning, got {:?}",
8851 result.warnings
8852 );
8853 }
8854
8855 #[test]
8856 fn count_free_params_includes_baseline_flags() {
8857 let data = u238_single_resonance();
8858 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
8859 let base = UnifiedFitConfig::new(
8860 energies,
8861 vec![data],
8862 vec!["U-238".into()],
8863 0.0,
8864 None,
8865 vec![0.001],
8866 )
8867 .unwrap();
8868 let n0 = count_free_params(&base);
8869
8870 let all_free = base
8871 .clone()
8872 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8873 assert_eq!(count_free_params(&all_free), n0 + 3);
8874
8875 let one_frozen = base
8876 .clone()
8877 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
8878 fit_b1: false,
8879 ..MultiplicativeBaselineConfig::default()
8880 });
8881 assert_eq!(count_free_params(&one_frozen), n0 + 2);
8882
8883 let frozen = base
8884 .clone()
8885 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
8886 fit_b0: false,
8887 fit_b1: false,
8888 fit_b2: false,
8889 ..MultiplicativeBaselineConfig::default()
8890 });
8891 assert_eq!(count_free_params(&frozen), n0);
8892 }
8893
8894 #[test]
8895 fn evaluate_jacobian_and_fisher_rejects_baseline() {
8896 let data = u238_single_resonance();
8897 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
8898 let n = energies.len();
8899 let config = UnifiedFitConfig::new(
8900 energies,
8901 vec![data],
8902 vec!["U-238".into()],
8903 0.0,
8904 None,
8905 vec![0.001],
8906 )
8907 .unwrap()
8908 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8909 let err = evaluate_jacobian_and_fisher(&config, &vec![1000.0; n], &vec![0.0; n])
8912 .err()
8913 .expect("research Fisher helper must reject a baseline config");
8914 assert!(err.to_string().contains("multiplicative"), "got: {err}");
8915 }
8916
8917 #[test]
8922 fn baseline_init_positivity_scoped_to_fit_window() {
8923 let data = u238_single_resonance();
8924 let energies: Vec<f64> = (0..61)
8928 .map(|i| 1e-3 * 10f64.powf(i as f64 / 10.0))
8929 .collect();
8930 let bl = MultiplicativeBaselineConfig {
8931 b0_init: 0.9,
8932 b2_init: -0.05,
8933 ..MultiplicativeBaselineConfig::default()
8934 };
8935 let base = UnifiedFitConfig::new(
8936 energies,
8937 vec![data],
8938 vec!["U-238".into()],
8939 0.0,
8940 None,
8941 vec![0.001],
8942 )
8943 .unwrap()
8944 .with_multiplicative_baseline(bl);
8945
8946 let err = validate_multiplicative_baseline(&base)
8948 .expect_err("negative B_init at unmasked edge bins must reject");
8949 assert!(err.to_string().contains("not strictly positive"), "{err}");
8950
8951 let windowed = base
8953 .with_fit_energy_range(Some((0.5, 2.0)))
8954 .expect("valid range");
8955 validate_multiplicative_baseline(&windowed)
8956 .expect("inits positive inside the fit window must be accepted");
8957 }
8958
8959 #[test]
8972 fn build_transmission_model_no_temp_matches_forward_model_bit_exact() {
8973 use nereids_physics::transmission::SampleParams;
8974 let data = u238_single_resonance();
8975 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8976 let density = 0.002;
8977
8978 let config = UnifiedFitConfig::new(
8980 energies.clone(),
8981 vec![data.clone()],
8982 vec!["U-238".into()],
8983 293.6,
8984 None,
8985 vec![density],
8986 )
8987 .unwrap();
8988 let model = build_transmission_model(&config, 1, None).unwrap();
8989 let t_model = model.evaluate(&[density]).unwrap();
8990 let sample = SampleParams::new(293.6, vec![(data.clone(), density)]).unwrap();
8991 let t_fwd = phys_transmission::forward_model(&energies, &sample, None).unwrap();
8992 assert_eq!(t_model.len(), t_fwd.len());
8993 for (i, (a, b)) in t_model.iter().zip(t_fwd.iter()).enumerate() {
8994 assert_eq!(
8995 a.to_bits(),
8996 b.to_bits(),
8997 "no-resolution arm, bin {i}: {a:e} != {b:e}"
8998 );
8999 }
9000
9001 let res = nereids_physics::resolution::ResolutionFunction::Gaussian(
9003 nereids_physics::resolution::ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap(),
9004 );
9005 let config_res = UnifiedFitConfig::new(
9006 energies.clone(),
9007 vec![data.clone()],
9008 vec!["U-238".into()],
9009 293.6,
9010 Some(res.clone()),
9011 vec![density],
9012 )
9013 .unwrap();
9014 let model_res = build_transmission_model(&config_res, 1, None).unwrap();
9015 let t_model_res = model_res.evaluate(&[density]).unwrap();
9016 let inst = InstrumentParams { resolution: res };
9017 let sample_res = SampleParams::new(293.6, vec![(data.clone(), density)]).unwrap();
9018 let t_fwd_res =
9019 phys_transmission::forward_model(&energies, &sample_res, Some(&inst)).unwrap();
9020 assert_eq!(t_model_res.len(), t_fwd_res.len());
9021 for (i, (a, b)) in t_model_res.iter().zip(t_fwd_res.iter()).enumerate() {
9022 assert_eq!(
9023 a.to_bits(),
9024 b.to_bits(),
9025 "Gaussian-resolution arm, bin {i}: {a:e} != {b:e}"
9026 );
9027 }
9028
9029 let working = phys_transmission::broadened_cross_sections_on_working_grid(
9030 &energies,
9031 std::slice::from_ref(&data),
9032 293.6,
9033 Some(&inst),
9034 None,
9035 )
9036 .unwrap();
9037 assert!(!working.layout.is_identity());
9038 let config_given = config_res
9039 .clone()
9040 .with_precomputed_cross_sections(PrecomputedXs::from(working));
9041 let t_given = build_transmission_model(&config_given, 1, None)
9042 .unwrap()
9043 .evaluate(&[density])
9044 .unwrap();
9045 for (i, (a, b)) in t_given.iter().zip(t_fwd_res.iter()).enumerate() {
9046 assert_eq!(
9047 a.to_bits(),
9048 b.to_bits(),
9049 "Gaussian-resolution arm with a supplied table, bin {i}: {a:e} != {b:e}"
9050 );
9051 }
9052
9053 let tab_text = "header\n---\n\
9054 5.0 0.0\n\
9055 -0.01 0.0\n\
9056 -0.005 0.5\n\
9057 0.0 1.0\n\
9058 0.005 0.5\n\
9059 0.01 0.0\n\
9060 \n\
9061 200.0 0.0\n\
9062 -0.02 0.0\n\
9063 -0.01 0.5\n\
9064 0.0 1.0\n\
9065 0.01 0.5\n\
9066 0.02 0.0\n";
9067 let tab = nereids_physics::resolution::TabulatedResolution::from_text(tab_text, 25.0)
9068 .expect("synthetic tabulated kernel parses");
9069 let res_tab = nereids_physics::resolution::ResolutionFunction::Tabulated(Arc::new(tab));
9070 let config_tab = UnifiedFitConfig::new(
9071 energies.clone(),
9072 vec![data.clone()],
9073 vec!["U-238".into()],
9074 293.6,
9075 Some(res_tab.clone()),
9076 vec![density],
9077 )
9078 .unwrap();
9079 let model_tab = build_transmission_model(&config_tab, 1, None).unwrap();
9080 let t_model_tab = model_tab.evaluate(&[density]).unwrap();
9081 let inst_tab = InstrumentParams {
9082 resolution: res_tab,
9083 };
9084 let sample_tab = SampleParams::new(293.6, vec![(data, density)]).unwrap();
9085 let t_fwd_tab =
9086 phys_transmission::forward_model(&energies, &sample_tab, Some(&inst_tab)).unwrap();
9087 assert_eq!(t_model_tab.len(), t_fwd_tab.len());
9088 for (i, (a, b)) in t_model_tab.iter().zip(t_fwd_tab.iter()).enumerate() {
9089 assert_eq!(
9090 a.to_bits(),
9091 b.to_bits(),
9092 "tabulated-resolution arm, bin {i}: {a:e} != {b:e}"
9093 );
9094 }
9095 let working_tab = phys_transmission::broadened_cross_sections_on_working_grid(
9096 &energies,
9097 std::slice::from_ref(&sample_tab.isotopes()[0].0),
9098 293.6,
9099 Some(&inst_tab),
9100 None,
9101 )
9102 .unwrap();
9103 assert!(!working_tab.layout.is_identity());
9104 let config_given_tab = config_tab
9105 .clone()
9106 .with_precomputed_cross_sections(PrecomputedXs::from(working_tab));
9107 let t_given_tab = build_transmission_model(&config_given_tab, 1, None)
9108 .unwrap()
9109 .evaluate(&[density])
9110 .unwrap();
9111 for (i, (a, b)) in t_given_tab.iter().zip(t_fwd_tab.iter()).enumerate() {
9112 assert_eq!(
9113 a.to_bits(),
9114 b.to_bits(),
9115 "tabulated-resolution arm with a supplied table, bin {i}: {a:e} != {b:e}"
9116 );
9117 }
9118 }
9119}