1use std::fmt;
11use std::sync::Arc;
12
13use nereids_core::constants::{EV_TO_JOULES, NEUTRON_MASS_KG};
14use nereids_endf::resonance::ResonanceData;
15use nereids_fitting::joint_poisson::{self, JointPoissonFitConfig, JointPoissonObjective};
16use nereids_fitting::lm::{self, FitModel, LmConfig, LmResult};
17use nereids_fitting::parameters::{FitParameter, ParameterSet};
18use nereids_fitting::poisson::{self, PoissonConfig};
19use nereids_fitting::transmission_model::{
20 EnergyScaleTransmissionModel, MultiplicativeBaselineModel, NormalizedTransmissionModel,
21 PrecomputedTransmissionModel, TransmissionFitModel,
22};
23use nereids_physics::resolution::ResolutionFunction;
24use nereids_physics::transmission::InstrumentParams;
25
26use crate::error::PipelineError;
27
28type PrecomputedWorkXs = (
33 Arc<Vec<Vec<f64>>>,
34 Arc<nereids_physics::transmission::WorkingGridLayout>,
35);
36
37#[derive(Debug, Clone)]
51pub struct BackgroundConfig {
52 pub anorm_init: f64,
54 pub back_a_init: f64,
56 pub back_b_init: f64,
58 pub back_c_init: f64,
60 pub back_d_init: f64,
65 pub back_f_init: f64,
70 pub fit_anorm: bool,
72 pub fit_back_a: bool,
74 pub fit_back_b: bool,
76 pub fit_back_c: bool,
78 pub fit_back_d: bool,
80 pub fit_back_f: bool,
82}
83
84impl Default for BackgroundConfig {
85 fn default() -> Self {
86 Self {
87 anorm_init: 1.0,
88 back_a_init: 0.0,
89 back_b_init: 0.0,
90 back_c_init: 0.0,
91 back_d_init: 0.01,
92 back_f_init: 1.0,
93 fit_anorm: true,
94 fit_back_a: true,
95 fit_back_b: true,
96 fit_back_c: true,
97 fit_back_d: false,
98 fit_back_f: false,
99 }
100 }
101}
102
103#[derive(Debug, Clone, Copy)]
108struct BackgroundIndices {
109 anorm: usize,
110 back_a: usize,
111 back_b: usize,
112 back_c: usize,
113 back_d: Option<usize>,
114 back_f: Option<usize>,
115}
116
117#[derive(Debug, Clone)]
149pub struct MultiplicativeBaselineConfig {
150 pub b0_init: f64,
152 pub b1_init: f64,
154 pub b2_init: f64,
156 pub fit_b0: bool,
159 pub fit_b1: bool,
161 pub fit_b2: bool,
163 pub spatial_global: bool,
169 pub b0_bounds: (f64, f64),
171 pub b1_bounds: (f64, f64),
173 pub b2_bounds: (f64, f64),
175}
176
177impl Default for MultiplicativeBaselineConfig {
178 fn default() -> Self {
179 Self {
180 b0_init: 1.0,
181 b1_init: 0.0,
182 b2_init: 0.0,
183 fit_b0: true,
184 fit_b1: true,
185 fit_b2: true,
186 spatial_global: true,
187 b0_bounds: (0.9, 1.1),
188 b1_bounds: (-0.05, 0.05),
189 b2_bounds: (-0.05, 0.05),
190 }
191 }
192}
193
194#[derive(Debug, Clone, Copy)]
197struct BaselineIndices {
198 b0: usize,
199 b1: usize,
200 b2: usize,
201}
202
203#[derive(Debug, Clone)]
215pub enum InputData {
216 Transmission {
220 transmission: Vec<f64>,
222 uncertainty: Vec<f64>,
224 },
225 Counts {
231 sample_counts: Vec<f64>,
233 open_beam_counts: Vec<f64>,
235 },
236 CountsWithNuisance {
240 sample_counts: Vec<f64>,
242 flux: Vec<f64>,
244 background: Vec<f64>,
246 },
247}
248
249impl InputData {
250 pub fn n_energies(&self) -> usize {
252 match self {
253 Self::Transmission { transmission, .. } => transmission.len(),
254 Self::Counts { sample_counts, .. } => sample_counts.len(),
255 Self::CountsWithNuisance { sample_counts, .. } => sample_counts.len(),
256 }
257 }
258
259 pub fn is_counts(&self) -> bool {
261 matches!(self, Self::Counts { .. } | Self::CountsWithNuisance { .. })
262 }
263}
264
265#[derive(Debug, Clone, Default)]
270pub enum SolverConfig {
271 LevenbergMarquardt(LmConfig),
273 PoissonKL(PoissonConfig),
288 #[default]
291 Auto,
292}
293
294#[derive(Debug, Clone)]
323pub struct CountsBackgroundConfig {
324 pub alpha_1_init: f64,
330 pub alpha_2_init: f64,
333 pub fit_alpha_1: bool,
337 pub fit_alpha_2: bool,
339 pub c: f64,
348}
349
350impl Default for CountsBackgroundConfig {
351 fn default() -> Self {
352 Self {
353 alpha_1_init: 1.0,
354 alpha_2_init: 1.0,
355 fit_alpha_1: false,
356 fit_alpha_2: false,
357 c: 1.0,
358 }
359 }
360}
361
362#[derive(Debug, Clone)]
369pub struct UnifiedFitConfig {
370 energies: Vec<f64>,
372 resonance_data: Vec<ResonanceData>,
373 isotope_names: Vec<String>,
374 temperature_k: f64,
375 resolution: Option<ResolutionFunction>,
376 initial_densities: Vec<f64>,
377 fit_temperature: bool,
378 compute_covariance: bool,
379 scale_by_chi2: bool,
385
386 solver: SolverConfig,
388
389 transmission_background: Option<BackgroundConfig>,
392 multiplicative_baseline: Option<MultiplicativeBaselineConfig>,
394 counts_background: Option<CountsBackgroundConfig>,
396
397 counts_enable_polish: Option<bool>,
404
405 precomputed_cross_sections: Option<Arc<Vec<Vec<f64>>>>,
407 precomputed_work_cross_sections: Option<PrecomputedWorkXs>,
425 precomputed_base_xs: Option<Arc<Vec<Vec<f64>>>>,
426 precomputed_resolution_plan: Option<Arc<nereids_physics::resolution::ResolutionPlan>>,
435 precomputed_sparse_cubature_plan:
446 Option<Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>>,
447 precomputed_sparse_scalar_plan: Option<Arc<nereids_physics::surrogate::ScalarSurrogatePlan>>,
451
452 fit_energy_scale: bool,
457 t0_init_us: f64,
459 l_scale_init: f64,
461 flight_path_m: f64,
463 tzero_jacobian_method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,
469 energy_scale_seed_enabled: bool,
477
478 fit_energy_range: Option<(f64, f64)>,
498
499 pub(crate) density_indices: Option<Vec<usize>>,
503 pub(crate) density_ratios: Option<Vec<f64>>,
506 n_density_params: Option<usize>,
509
510 density_free: Option<Vec<bool>>,
527}
528
529impl UnifiedFitConfig {
530 pub fn new(
532 energies: Vec<f64>,
533 resonance_data: Vec<ResonanceData>,
534 isotope_names: Vec<String>,
535 temperature_k: f64,
536 resolution: Option<ResolutionFunction>,
537 initial_densities: Vec<f64>,
538 ) -> Result<Self, FitConfigError> {
539 if energies.is_empty() {
540 return Err(FitConfigError::EmptyEnergies);
541 }
542 if resonance_data.is_empty() {
543 return Err(FitConfigError::EmptyResonanceData);
544 }
545 if initial_densities.len() != resonance_data.len() {
546 return Err(FitConfigError::DensityCountMismatch {
547 densities: initial_densities.len(),
548 isotopes: resonance_data.len(),
549 });
550 }
551 if isotope_names.len() != resonance_data.len() {
552 return Err(FitConfigError::NameCountMismatch {
553 names: isotope_names.len(),
554 isotopes: resonance_data.len(),
555 });
556 }
557 if !temperature_k.is_finite() {
558 return Err(FitConfigError::NonFiniteTemperature(temperature_k));
559 }
560 if temperature_k < 0.0 {
561 return Err(FitConfigError::NegativeTemperature(temperature_k));
562 }
563 Ok(Self {
564 energies,
565 resonance_data,
566 isotope_names,
567 temperature_k,
568 resolution,
569 initial_densities,
570 fit_temperature: false,
571 compute_covariance: true,
572 scale_by_chi2: false,
573 solver: SolverConfig::Auto,
574 transmission_background: None,
575 multiplicative_baseline: None,
576 counts_background: None,
577 counts_enable_polish: None,
578 precomputed_cross_sections: None,
579 precomputed_work_cross_sections: None,
580 precomputed_base_xs: None,
581 precomputed_resolution_plan: None,
582 precomputed_sparse_cubature_plan: None,
583 precomputed_sparse_scalar_plan: None,
584 fit_energy_scale: false,
585 t0_init_us: 0.0,
586 l_scale_init: 1.0,
587 flight_path_m: 25.0,
588 density_indices: None,
589 density_ratios: None,
590 n_density_params: None,
591 density_free: None,
592 tzero_jacobian_method: None,
593 energy_scale_seed_enabled: true,
594 fit_energy_range: None,
595 })
596 }
597
598 #[must_use]
604 pub fn with_tzero_jacobian_method(
605 mut self,
606 method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,
607 ) -> Self {
608 self.tzero_jacobian_method = method;
609 self
610 }
611
612 #[must_use]
619 pub fn with_energy_scale_seed(mut self, enabled: bool) -> Self {
620 self.energy_scale_seed_enabled = enabled;
621 self
622 }
623
624 #[must_use]
627 pub fn with_solver(mut self, solver: SolverConfig) -> Self {
628 self.solver = solver;
629 self
630 }
631
632 #[must_use]
633 pub fn with_fit_temperature(mut self, v: bool) -> Self {
634 self.fit_temperature = v;
635 self
636 }
637
638 #[must_use]
639 pub fn with_compute_covariance(mut self, v: bool) -> Self {
640 self.compute_covariance = v;
641 self
642 }
643
644 #[must_use]
664 pub fn with_scale_by_chi2(mut self, v: bool) -> Self {
665 self.scale_by_chi2 = v;
666 self
667 }
668
669 #[must_use]
675 pub fn with_energy_scale(
676 mut self,
677 t0_init_us: f64,
678 l_scale_init: f64,
679 flight_path_m: f64,
680 ) -> Self {
681 self.fit_energy_scale = true;
682 self.t0_init_us = t0_init_us;
683 self.l_scale_init = l_scale_init;
684 self.flight_path_m = flight_path_m;
685 self
686 }
687
688 pub fn with_fit_energy_range(
701 mut self,
702 range: Option<(f64, f64)>,
703 ) -> Result<Self, FitConfigError> {
704 if let Some((lo, hi)) = range {
705 if !lo.is_finite() || !hi.is_finite() {
706 return Err(FitConfigError::InvalidFitEnergyRange(
707 "fit_energy_range bounds must be finite",
708 ));
709 }
710 if lo >= hi {
711 return Err(FitConfigError::InvalidFitEnergyRange(
712 "fit_energy_range min must be strictly less than max",
713 ));
714 }
715 }
716 self.fit_energy_range = range;
717 Ok(self)
718 }
719
720 #[must_use]
721 pub fn with_transmission_background(mut self, bg: BackgroundConfig) -> Self {
722 self.transmission_background = Some(bg);
723 self
724 }
725
726 #[must_use]
733 pub fn with_multiplicative_baseline(mut self, bl: MultiplicativeBaselineConfig) -> Self {
734 self.multiplicative_baseline = Some(bl);
735 self
736 }
737
738 #[must_use]
739 pub fn with_counts_background(mut self, bg: CountsBackgroundConfig) -> Self {
740 self.counts_background = Some(bg);
741 self
742 }
743
744 #[must_use]
749 pub fn with_counts_enable_polish(mut self, v: Option<bool>) -> Self {
750 self.counts_enable_polish = v;
751 self
752 }
753
754 #[must_use]
755 pub fn with_precomputed_cross_sections(mut self, xs: Arc<Vec<Vec<f64>>>) -> Self {
756 self.precomputed_cross_sections = Some(xs);
757 self.precomputed_sparse_cubature_plan = None;
762 self.precomputed_sparse_scalar_plan = None;
763 self.precomputed_work_cross_sections = None;
767 self
768 }
769
770 #[must_use]
789 pub fn with_precomputed_work_cross_sections(
790 mut self,
791 xs: Arc<Vec<Vec<f64>>>,
792 layout: Arc<nereids_physics::transmission::WorkingGridLayout>,
793 ) -> Self {
794 self.precomputed_work_cross_sections = Some((xs, layout));
795 self
796 }
797
798 #[must_use]
799 pub fn with_precomputed_base_xs(mut self, xs: Arc<Vec<Vec<f64>>>) -> Self {
800 self.precomputed_base_xs = Some(xs);
801 self.precomputed_sparse_cubature_plan = None;
805 self.precomputed_sparse_scalar_plan = None;
806 self
807 }
808
809 #[must_use]
817 pub fn with_precomputed_resolution_plan(
818 mut self,
819 plan: Arc<nereids_physics::resolution::ResolutionPlan>,
820 ) -> Self {
821 self.precomputed_resolution_plan = Some(plan);
822 self
823 }
824
825 #[must_use]
843 pub fn with_precomputed_sparse_cubature_plan(
844 mut self,
845 plan: Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>,
846 ) -> Self {
847 self.precomputed_sparse_cubature_plan = Some(plan);
848 self
849 }
850
851 #[must_use]
857 pub fn with_precomputed_sparse_scalar_plan(
858 mut self,
859 plan: Arc<nereids_physics::surrogate::ScalarSurrogatePlan>,
860 ) -> Self {
861 self.precomputed_sparse_scalar_plan = Some(plan);
862 self
863 }
864
865 pub fn with_groups(
882 mut self,
883 groups: &[(&nereids_core::types::IsotopeGroup, &[ResonanceData])],
884 initial_densities: Vec<f64>,
885 ) -> Result<Self, FitConfigError> {
886 if self.density_free.is_some() {
887 return Err(FitConfigError::DensityFreezeBeforeGroups);
888 }
889 if groups.is_empty() {
890 return Err(FitConfigError::EmptyResonanceData);
891 }
892 if initial_densities.len() != groups.len() {
893 return Err(FitConfigError::DensityCountMismatch {
894 densities: initial_densities.len(),
895 isotopes: groups.len(),
896 });
897 }
898 let mut all_resonance_data = Vec::new();
899 let mut all_indices = Vec::new();
900 let mut all_ratios = Vec::new();
901 let mut names = Vec::new();
902 for (g_idx, (group, rd_list)) in groups.iter().enumerate() {
903 if rd_list.len() != group.n_members() {
904 return Err(FitConfigError::GroupMemberCountMismatch {
905 group_name: group.name().to_string(),
906 rd_count: rd_list.len(),
907 member_count: group.n_members(),
908 });
909 }
910 names.push(group.name().to_string());
911 for ((isotope, ratio), rd) in group.members().iter().zip(rd_list.iter()) {
912 if rd.isotope != *isotope {
914 return Err(FitConfigError::GroupMemberIsotopeMismatch {
915 group_name: group.name().to_string(),
916 expected_z: isotope.z(),
917 expected_a: isotope.a(),
918 got_z: rd.isotope.z(),
919 got_a: rd.isotope.a(),
920 });
921 }
922 all_resonance_data.push(rd.clone());
923 all_indices.push(g_idx);
924 all_ratios.push(*ratio);
925 }
926 }
927 self.resonance_data = all_resonance_data;
928 self.isotope_names = names;
929 self.initial_densities = initial_densities;
930 self.n_density_params = Some(groups.len());
931 self.density_indices = Some(all_indices);
932 self.density_ratios = Some(all_ratios);
933 self.precomputed_cross_sections = None;
939 self.precomputed_work_cross_sections = None;
940 self.precomputed_base_xs = None;
941 self.precomputed_sparse_cubature_plan = None;
947 self.precomputed_sparse_scalar_plan = None;
948 Ok(self)
949 }
950
951 pub fn precomputed_sparse_cubature_plan(
958 &self,
959 ) -> Option<&Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>> {
960 self.precomputed_sparse_cubature_plan.as_ref()
961 }
962
963 pub fn precomputed_sparse_scalar_plan(
965 &self,
966 ) -> Option<&Arc<nereids_physics::surrogate::ScalarSurrogatePlan>> {
967 self.precomputed_sparse_scalar_plan.as_ref()
968 }
969
970 pub fn energies(&self) -> &[f64] {
971 &self.energies
972 }
973 pub fn resonance_data(&self) -> &[ResonanceData] {
974 &self.resonance_data
975 }
976 pub fn isotope_names(&self) -> &[String] {
977 &self.isotope_names
978 }
979 pub fn temperature_k(&self) -> f64 {
980 self.temperature_k
981 }
982 pub fn resolution(&self) -> Option<&ResolutionFunction> {
983 self.resolution.as_ref()
984 }
985 pub fn initial_densities(&self) -> &[f64] {
986 &self.initial_densities
987 }
988 pub fn solver(&self) -> &SolverConfig {
989 &self.solver
990 }
991 pub fn fit_temperature(&self) -> bool {
992 self.fit_temperature
993 }
994 pub fn transmission_background(&self) -> Option<&BackgroundConfig> {
995 self.transmission_background.as_ref()
996 }
997 pub fn multiplicative_baseline(&self) -> Option<&MultiplicativeBaselineConfig> {
999 self.multiplicative_baseline.as_ref()
1000 }
1001 pub fn counts_background(&self) -> Option<&CountsBackgroundConfig> {
1002 self.counts_background.as_ref()
1003 }
1004 pub fn counts_enable_polish(&self) -> Option<bool> {
1006 self.counts_enable_polish
1007 }
1008 pub fn fit_energy_scale(&self) -> bool {
1011 self.fit_energy_scale
1012 }
1013 pub fn scale_by_chi2(&self) -> bool {
1016 self.scale_by_chi2
1017 }
1018 pub fn flight_path_m(&self) -> f64 {
1020 self.flight_path_m
1021 }
1022 pub fn fit_energy_range(&self) -> Option<(f64, f64)> {
1026 self.fit_energy_range
1027 }
1028 pub fn baseline_reference_energy(&self) -> f64 {
1036 let mask = nereids_fitting::active_mask::build_active_mask(
1037 self.energies(),
1038 self.fit_energy_range(),
1039 );
1040 nereids_fitting::transmission_model::baseline_reference_energy_active(
1041 self.energies(),
1042 mask.as_deref(),
1043 )
1044 }
1045 pub fn precomputed_cross_sections(&self) -> Option<&Arc<Vec<Vec<f64>>>> {
1046 self.precomputed_cross_sections.as_ref()
1047 }
1048 pub fn n_density_params(&self) -> usize {
1050 self.n_density_params.unwrap_or(self.resonance_data.len())
1051 }
1052
1053 #[must_use]
1066 pub fn with_fix_densities(mut self, fix: bool) -> Self {
1067 self.density_free = if fix {
1068 Some(vec![false; self.n_density_params()])
1069 } else {
1070 None
1071 };
1072 self
1073 }
1074
1075 pub fn with_density_free(mut self, free: Vec<bool>) -> Result<Self, FitConfigError> {
1089 let n = self.n_density_params();
1090 if free.len() != n {
1091 return Err(FitConfigError::DensityCountMismatch {
1092 densities: free.len(),
1093 isotopes: n,
1094 });
1095 }
1096 self.density_free = if free.iter().all(|&f| f) {
1099 None
1100 } else {
1101 Some(free)
1102 };
1103 Ok(self)
1104 }
1105
1106 fn density_is_fixed(&self, i: usize) -> bool {
1109 self.density_free
1110 .as_ref()
1111 .is_some_and(|mask| !mask.get(i).copied().unwrap_or(true))
1112 }
1113
1114 fn n_free_density_params(&self) -> usize {
1117 (0..self.n_density_params())
1118 .filter(|&i| !self.density_is_fixed(i))
1119 .count()
1120 }
1121
1122 pub(crate) fn effective_solver(&self, input: &InputData) -> SolverConfig {
1124 match &self.solver {
1125 SolverConfig::Auto => {
1126 if input.is_counts() {
1127 SolverConfig::PoissonKL(PoissonConfig::default())
1128 } else {
1129 SolverConfig::LevenbergMarquardt(LmConfig::default())
1130 }
1131 }
1132 other => other.clone(),
1133 }
1134 }
1135}
1136
1137pub fn fit_spectrum_typed(
1150 input: &InputData,
1151 config: &UnifiedFitConfig,
1152) -> Result<SpectrumFitResult, PipelineError> {
1153 let n_e = config.energies().len();
1154
1155 if config.fit_temperature && config.temperature_k < 1.0 {
1157 return Err(PipelineError::InvalidParameter(format!(
1158 "temperature must be >= 1.0 K when fit_temperature is true, got {}",
1159 config.temperature_k,
1160 )));
1161 }
1162
1163 if count_free_params(config) == 0 {
1169 return Err(PipelineError::InvalidParameter(
1170 "no free parameters to fit: all densities are frozen and no other \
1171 parameter is free — free at least one density (with_density_free) \
1172 or enable fit_temperature / energy-scale / background"
1173 .to_string(),
1174 ));
1175 }
1176
1177 if input.n_energies() != n_e {
1179 return Err(PipelineError::ShapeMismatch(format!(
1180 "input data has {} energy bins but config.energies has {}",
1181 input.n_energies(),
1182 n_e,
1183 )));
1184 }
1185
1186 match input {
1188 InputData::Transmission {
1189 transmission,
1190 uncertainty,
1191 } => {
1192 if uncertainty.len() != transmission.len() {
1193 return Err(PipelineError::ShapeMismatch(format!(
1194 "uncertainty length {} != transmission length {}",
1195 uncertainty.len(),
1196 transmission.len(),
1197 )));
1198 }
1199 }
1200 InputData::Counts {
1201 sample_counts,
1202 open_beam_counts,
1203 } => {
1204 if open_beam_counts.len() != sample_counts.len() {
1205 return Err(PipelineError::ShapeMismatch(format!(
1206 "open_beam_counts length {} != sample_counts length {}",
1207 open_beam_counts.len(),
1208 sample_counts.len(),
1209 )));
1210 }
1211 }
1212 InputData::CountsWithNuisance {
1213 sample_counts,
1214 flux,
1215 background,
1216 } => {
1217 if flux.len() != sample_counts.len() {
1218 return Err(PipelineError::ShapeMismatch(format!(
1219 "flux length {} != sample_counts length {}",
1220 flux.len(),
1221 sample_counts.len(),
1222 )));
1223 }
1224 if background.len() != sample_counts.len() {
1225 return Err(PipelineError::ShapeMismatch(format!(
1226 "background length {} != sample_counts length {}",
1227 background.len(),
1228 sample_counts.len(),
1229 )));
1230 }
1231 }
1232 }
1233
1234 validate_precomputed_cross_sections(config)?;
1238
1239 let effective_solver = config.effective_solver(input);
1240
1241 match (input, &effective_solver) {
1242 (
1244 InputData::Transmission {
1245 transmission,
1246 uncertainty,
1247 },
1248 SolverConfig::LevenbergMarquardt(lm_cfg),
1249 ) => fit_transmission_lm(transmission, uncertainty, config, lm_cfg),
1250
1251 (
1253 InputData::Transmission {
1254 transmission,
1255 uncertainty,
1256 },
1257 SolverConfig::PoissonKL(poisson_cfg),
1258 ) => fit_transmission_poisson(transmission, uncertainty, config, poisson_cfg),
1259
1260 (
1270 InputData::Counts {
1271 sample_counts,
1272 open_beam_counts,
1273 },
1274 SolverConfig::PoissonKL(poisson_cfg),
1275 ) => {
1276 let bg = vec![0.0f64; n_e];
1277 fit_counts_joint_poisson(
1278 sample_counts,
1279 open_beam_counts,
1280 &bg,
1281 config,
1282 &poisson_to_joint_poisson_config(poisson_cfg, config),
1283 )
1284 }
1285
1286 (
1288 InputData::CountsWithNuisance {
1289 sample_counts,
1290 flux,
1291 background,
1292 },
1293 SolverConfig::PoissonKL(poisson_cfg),
1294 ) => fit_counts_joint_poisson(
1295 sample_counts,
1296 flux,
1297 background,
1298 config,
1299 &poisson_to_joint_poisson_config(poisson_cfg, config),
1300 ),
1301
1302 (
1311 InputData::Counts {
1312 sample_counts,
1313 open_beam_counts,
1314 },
1315 SolverConfig::LevenbergMarquardt(lm_cfg),
1316 ) => {
1317 let (transmission, uncertainty) =
1318 counts_to_transmission(sample_counts, open_beam_counts);
1319 fit_transmission_lm(&transmission, &uncertainty, config, lm_cfg)
1320 }
1321
1322 (InputData::CountsWithNuisance { .. }, SolverConfig::LevenbergMarquardt(_)) => {
1324 Err(PipelineError::InvalidParameter(
1325 "CountsWithNuisance requires a counts-domain solver (LM cannot use nuisance parameters)"
1326 .into(),
1327 ))
1328 }
1329
1330 (_, SolverConfig::Auto) => unreachable!("Auto should be resolved before dispatch"),
1332 }
1333}
1334
1335fn poisson_to_joint_poisson_config(
1341 poisson_cfg: &PoissonConfig,
1342 config: &UnifiedFitConfig,
1343) -> JointPoissonFitConfig {
1344 let mut jp_cfg = JointPoissonFitConfig {
1345 max_iter: poisson_cfg.max_iter,
1346 scale_by_chi2: config.scale_by_chi2,
1347 ..Default::default()
1348 };
1349 if let Some(override_val) = config.counts_enable_polish() {
1350 jp_cfg.enable_polish = override_val;
1351 }
1352 jp_cfg
1353}
1354
1355fn counts_to_transmission(sample: &[f64], open_beam: &[f64]) -> (Vec<f64>, Vec<f64>) {
1369 let transmission: Vec<f64> = sample
1370 .iter()
1371 .zip(open_beam.iter())
1372 .map(|(&s, &ob)| if ob > 0.0 { s / ob } else { 0.0 })
1373 .collect();
1374 let uncertainty: Vec<f64> = sample
1375 .iter()
1376 .zip(open_beam.iter())
1377 .map(|(&s, &ob)| {
1378 if ob <= 0.0 {
1379 1e30
1381 } else if s <= 0.0 {
1382 1e10
1384 } else {
1385 s.max(1.0).sqrt() / ob.max(1e-10)
1386 }
1387 })
1388 .collect();
1389 (transmission, uncertainty)
1390}
1391
1392fn fit_transmission_lm(
1394 measured_t: &[f64],
1395 sigma: &[f64],
1396 config: &UnifiedFitConfig,
1397 lm_config: &LmConfig,
1398) -> Result<SpectrumFitResult, PipelineError> {
1399 let n_density_params = config.n_density_params();
1400
1401 let mut param_vec = build_density_params(config);
1403
1404 let temperature_index = append_temperature_param(&mut param_vec, config);
1405 let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);
1406
1407 seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, measured_t, config);
1411
1412 if let Some(bg) = config.transmission_background.as_ref() {
1415 validate_transmission_background(bg)?;
1416 }
1417 let bg_indices = config
1418 .transmission_background
1419 .as_ref()
1420 .map(|bg| append_background_params(&mut param_vec, bg));
1421
1422 validate_multiplicative_baseline(config)?;
1425 let bl_indices = config
1426 .multiplicative_baseline
1427 .as_ref()
1428 .map(|bl| append_multiplicative_baseline_params(&mut param_vec, bl));
1429
1430 let mut params = ParameterSet::new(param_vec);
1431 let mut lm_cfg = lm_config.clone();
1432 lm_cfg.compute_covariance = config.compute_covariance;
1433
1434 let model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1436 build_energy_scale_transmission_model(config, t0_idx, ls_idx, temperature_index)?
1437 } else {
1438 build_transmission_model(config, n_density_params, temperature_index)?
1439 };
1440
1441 let active_mask = nereids_fitting::active_mask::build_active_mask(
1445 config.energies(),
1446 config.fit_energy_range(),
1447 );
1448
1449 if let Some((e_min, e_max)) = config.fit_energy_range() {
1460 let n_active = nereids_fitting::active_mask::active_count(
1461 active_mask.as_deref(),
1462 config.energies().len(),
1463 );
1464 let required = required_active_bins(config);
1465 if n_active < required {
1466 return Err(PipelineError::InvalidParameter(format!(
1467 "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
1468 on the configured energy grid; at least {required} active bin(s) are \
1469 required for LM transmission fitting with {n_free} free parameter(s) \
1470 (underdetermined when n_active < n_free)",
1471 n_free = count_free_params(config),
1472 )));
1473 }
1474 }
1475
1476 let mut stacked: Box<dyn FitModel> = model;
1482 if let Some(bi) = bg_indices {
1483 stacked = if let (Some(di), Some(fi)) = (bi.back_d, bi.back_f) {
1484 Box::new(NormalizedTransmissionModel::new_with_exponential(
1485 stacked,
1486 config.energies(),
1487 bi.anorm,
1488 bi.back_a,
1489 bi.back_b,
1490 bi.back_c,
1491 di,
1492 fi,
1493 ))
1494 } else {
1495 Box::new(NormalizedTransmissionModel::new(
1496 stacked,
1497 config.energies(),
1498 bi.anorm,
1499 bi.back_a,
1500 bi.back_b,
1501 bi.back_c,
1502 ))
1503 };
1504 }
1505 if let Some(bli) = bl_indices {
1506 stacked = Box::new(
1507 MultiplicativeBaselineModel::new(
1508 stacked,
1509 config.energies(),
1510 nereids_fitting::transmission_model::baseline_reference_energy_active(
1511 config.energies(),
1512 active_mask.as_deref(),
1513 ),
1514 bli.b0,
1515 bli.b1,
1516 bli.b2,
1517 )
1518 .with_active_mask(active_mask.as_deref()),
1520 );
1521 }
1522 let result = lm::levenberg_marquardt_with_mask(
1523 &*stacked,
1524 measured_t,
1525 sigma,
1526 &mut params,
1527 &lm_cfg,
1528 active_mask.as_deref(),
1529 )?;
1530
1531 let free_indices = params.free_indices();
1532 let mut sr = extract_result(
1533 config,
1534 &result,
1535 n_density_params,
1536 &free_indices,
1537 bg_indices,
1538 bl_indices,
1539 )?;
1540
1541 if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1543 sr.t0_us = Some(result.params[t0_idx]);
1544 sr.l_scale = Some(result.params[ls_idx]);
1545 sr.energy_scale_flight_path_m = Some(config.flight_path_m);
1546 }
1547
1548 Ok(sr)
1549}
1550
1551fn fit_transmission_poisson(
1558 measured_t: &[f64],
1559 sigma: &[f64],
1560 config: &UnifiedFitConfig,
1561 poisson_cfg: &PoissonConfig,
1562) -> Result<SpectrumFitResult, PipelineError> {
1563 if config.fit_energy_range().is_some() {
1571 return Err(PipelineError::InvalidParameter(
1572 "fit_energy_range is not supported for the transmission + \
1573 Poisson-KL solver path. Use joint-Poisson (provide sample + \
1574 open-beam counts) or switch to the LM transmission solver."
1575 .into(),
1576 ));
1577 }
1578
1579 let mut poisson_cfg = poisson_cfg.clone();
1580 poisson_cfg.compute_covariance = config.compute_covariance;
1581 let poisson_cfg = &poisson_cfg;
1582
1583 let n_density_params = config.n_density_params();
1584 let mut param_vec = build_density_params(config);
1585
1586 let temperature_index = append_temperature_param(&mut param_vec, config);
1587 let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);
1588
1589 seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, measured_t, config);
1591
1592 if let Some(bg) = config.transmission_background.as_ref() {
1595 validate_transmission_background(bg)?;
1596 }
1597 let bg_indices = config
1598 .transmission_background
1599 .as_ref()
1600 .map(|bg| append_background_params(&mut param_vec, bg));
1601
1602 validate_multiplicative_baseline(config)?;
1605 let bl_indices = config
1606 .multiplicative_baseline
1607 .as_ref()
1608 .map(|bl| append_multiplicative_baseline_params(&mut param_vec, bl));
1609
1610 let mut params = ParameterSet::new(param_vec);
1611
1612 let model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1614 build_energy_scale_transmission_model(config, t0_idx, ls_idx, temperature_index)?
1615 } else {
1616 build_transmission_model(config, n_density_params, temperature_index)?
1617 };
1618
1619 let mut stacked: Box<dyn FitModel> = model;
1623 if let Some(bi) = bg_indices {
1624 stacked = if let (Some(di), Some(fi)) = (bi.back_d, bi.back_f) {
1625 Box::new(NormalizedTransmissionModel::new_with_exponential(
1626 stacked,
1627 config.energies(),
1628 bi.anorm,
1629 bi.back_a,
1630 bi.back_b,
1631 bi.back_c,
1632 di,
1633 fi,
1634 ))
1635 } else {
1636 Box::new(NormalizedTransmissionModel::new(
1637 stacked,
1638 config.energies(),
1639 bi.anorm,
1640 bi.back_a,
1641 bi.back_b,
1642 bi.back_c,
1643 ))
1644 };
1645 }
1646 if let Some(bli) = bl_indices {
1647 stacked = Box::new(MultiplicativeBaselineModel::new(
1651 stacked,
1652 config.energies(),
1653 config.baseline_reference_energy(),
1654 bli.b0,
1655 bli.b1,
1656 bli.b2,
1657 ));
1658 }
1659 let pr = poisson::poisson_fit(&*stacked, measured_t, &mut params, poisson_cfg)?;
1660 let result = poisson_to_lm_result(
1661 &*stacked,
1662 measured_t,
1663 sigma,
1664 &pr,
1665 ¶ms,
1666 config.scale_by_chi2,
1667 )?;
1668
1669 let free_indices = params.free_indices();
1670 let mut sr = extract_result(
1671 config,
1672 &result,
1673 n_density_params,
1674 &free_indices,
1675 bg_indices,
1676 bl_indices,
1677 )?;
1678
1679 if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1687 sr.t0_us = Some(result.params[t0_idx]);
1688 sr.l_scale = Some(result.params[ls_idx]);
1689 sr.energy_scale_flight_path_m = Some(config.flight_path_m);
1690 }
1691
1692 Ok(sr)
1693}
1694
1695fn fit_counts_joint_poisson(
1714 sample_counts: &[f64],
1715 flux: &[f64],
1716 detector_background: &[f64],
1717 config: &UnifiedFitConfig,
1718 jp_cfg: &JointPoissonFitConfig,
1719) -> Result<SpectrumFitResult, PipelineError> {
1720 if let Some(bg) = config.counts_background()
1722 && (bg.fit_alpha_1 || bg.fit_alpha_2)
1723 {
1724 return Err(PipelineError::InvalidParameter(
1725 "joint-Poisson solver does not support fit_alpha_1/fit_alpha_2: \
1726 the profile lambda-hat absorbs the global flux scale (alpha_1 redundant); \
1727 alpha_2 / B_det wiring is not yet implemented."
1728 .into(),
1729 ));
1730 }
1731 if detector_background.iter().any(|&v| v.abs() > 1e-12) {
1732 return Err(PipelineError::InvalidParameter(
1733 "joint-Poisson solver with non-zero detector_background is not yet supported \
1734 (B_det wiring is deferred)."
1735 .into(),
1736 ));
1737 }
1738
1739 if let Some(bg) = config.transmission_background.as_ref() {
1741 if bg.fit_back_d || bg.fit_back_f {
1742 return Err(PipelineError::InvalidParameter(
1743 "joint-Poisson solver does not support the BackD/BackF exponential \
1744 tail (support is deferred)."
1745 .into(),
1746 ));
1747 }
1748 if (bg.fit_back_b || bg.fit_back_c) && !bg.fit_back_a {
1749 return Err(PipelineError::InvalidParameter(
1750 "joint-Poisson transmission_background: B_A (fit_back_a) must be \
1751 enabled whenever any of B_B / B_C is enabled (A_n alone cannot \
1752 absorb a constant offset — benchmarked at −23% density bias)."
1753 .into(),
1754 ));
1755 }
1756 }
1757
1758 let c = config.counts_background().map(|b| b.c).unwrap_or(1.0);
1759 if !(c.is_finite() && c > 0.0) {
1760 return Err(PipelineError::InvalidParameter(format!(
1761 "joint-Poisson solver requires finite c > 0 in CountsBackgroundConfig, got {c}",
1762 )));
1763 }
1764
1765 let n_density_params = config.n_density_params();
1770 let mut param_vec = build_density_params(config);
1771
1772 let temperature_index = append_temperature_param(&mut param_vec, config);
1773 let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);
1774
1775 if energy_scale_indices.is_some() {
1780 let t_proxy: Vec<f64> = sample_counts
1781 .iter()
1782 .zip(flux.iter())
1783 .zip(detector_background.iter())
1784 .map(|((&s, &f), &b)| {
1785 let den = f - b;
1786 if den > 0.0 {
1787 ((s - b).max(0.0) / den).min(2.0)
1788 } else {
1789 1.0
1790 }
1791 })
1792 .collect();
1793 seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, &t_proxy, config);
1794 }
1795
1796 let bg_indices = config
1800 .transmission_background
1801 .as_ref()
1802 .map(|bg| append_background_params(&mut param_vec, bg));
1803
1804 validate_multiplicative_baseline(config)?;
1809 let bl_indices = config
1810 .multiplicative_baseline
1811 .as_ref()
1812 .map(|bl| append_multiplicative_baseline_params(&mut param_vec, bl));
1813
1814 let mut params = ParameterSet::new(param_vec);
1815
1816 let t_model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
1818 build_energy_scale_transmission_model(config, t0_idx, ls_idx, temperature_index)?
1819 } else {
1820 build_transmission_model(config, n_density_params, temperature_index)?
1821 };
1822
1823 let active_mask = nereids_fitting::active_mask::build_active_mask(
1834 config.energies(),
1835 config.fit_energy_range(),
1836 );
1837 let active_mask_slice = active_mask.as_deref();
1838
1839 if let Some((e_min, e_max)) = config.fit_energy_range() {
1849 let n_active =
1850 nereids_fitting::active_mask::active_count(active_mask_slice, config.energies().len());
1851 let required = required_active_bins(config);
1852 if n_active < required {
1853 return Err(PipelineError::InvalidParameter(format!(
1854 "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
1855 on the configured energy grid; at least {required} active bin(s) are \
1856 required for joint-Poisson fitting with {n_free} free parameter(s) \
1857 (underdetermined when n_active < n_free)",
1858 n_free = count_free_params(config),
1859 )));
1860 }
1861 }
1862
1863 let mut stacked: Box<dyn FitModel> = t_model;
1867 if let Some(bi) = bg_indices {
1868 stacked = Box::new(NormalizedTransmissionModel::new(
1869 stacked,
1870 config.energies(),
1871 bi.anorm,
1872 bi.back_a,
1873 bi.back_b,
1874 bi.back_c,
1875 ));
1876 }
1877 if let Some(bli) = bl_indices {
1878 stacked = Box::new(
1879 MultiplicativeBaselineModel::new(
1880 stacked,
1881 config.energies(),
1882 nereids_fitting::transmission_model::baseline_reference_energy_active(
1883 config.energies(),
1884 active_mask_slice,
1885 ),
1886 bli.b0,
1887 bli.b1,
1888 bli.b2,
1889 )
1890 .with_active_mask(active_mask_slice),
1892 );
1893 }
1894 let objective = JointPoissonObjective {
1895 model: &*stacked,
1896 o: flux,
1897 s: sample_counts,
1898 c,
1899 active_mask: active_mask_slice,
1900 };
1901 let mut cfg = jp_cfg.clone();
1902 cfg.compute_covariance = config.compute_covariance;
1903 cfg.scale_by_chi2 = config.scale_by_chi2;
1904 let result = joint_poisson::joint_poisson_fit(&objective, &mut params, &cfg)
1912 .map_err(PipelineError::Fitting)?;
1913
1914 let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();
1916
1917 let (uncertainties, temperature_k_unc) = if let Some(ref unc_all) = result.uncertainties {
1918 let free_idx = params.free_indices();
1922 let dens_unc: Vec<f64> = (0..n_density_params)
1923 .map(|i| free_uncertainty(&free_idx, unc_all, i).unwrap_or(f64::NAN))
1924 .collect();
1925 let t_unc = temperature_index.and_then(|idx| free_uncertainty(&free_idx, unc_all, idx));
1926 (Some(dens_unc), t_unc)
1927 } else {
1928 (None, None)
1929 };
1930 let fitted_temp = temperature_index.map(|idx| result.params[idx]);
1931
1932 let converged = result.gn_converged || result.polish_converged;
1937
1938 let (anorm_out, bg_abc_out) = if let Some(bi) = bg_indices {
1943 (
1944 result.params[bi.anorm],
1945 [
1946 result.params[bi.back_a],
1947 result.params[bi.back_b],
1948 result.params[bi.back_c],
1949 ],
1950 )
1951 } else {
1952 (1.0, [0.0, 0.0, 0.0])
1953 };
1954
1955 let (baseline_out, baseline_e_ref_out) = if let Some(bli) = bl_indices {
1957 (
1958 Some([
1959 result.params[bli.b0],
1960 result.params[bli.b1],
1961 result.params[bli.b2],
1962 ]),
1963 Some(
1964 nereids_fitting::transmission_model::baseline_reference_energy_active(
1965 config.energies(),
1966 active_mask_slice,
1967 ),
1968 ),
1969 )
1970 } else {
1971 (None, None)
1972 };
1973
1974 Ok(SpectrumFitResult {
1975 densities,
1976 uncertainties,
1977 reduced_chi_squared: result.deviance_per_dof,
1981 converged,
1982 iterations: result.gn_iterations + result.polish_iterations,
1983 temperature_k: fitted_temp,
1984 temperature_k_unc,
1985 anorm: anorm_out,
1986 background: bg_abc_out,
1987 back_d: None,
1992 back_f: None,
1993 t0_us: energy_scale_indices.map(|(t0_idx, _)| result.params[t0_idx]),
1994 l_scale: energy_scale_indices.map(|(_, ls_idx)| result.params[ls_idx]),
1995 energy_scale_flight_path_m: energy_scale_indices.map(|_| config.flight_path_m),
1996 deviance_per_dof: Some(result.deviance_per_dof),
1997 baseline: baseline_out,
1998 baseline_e_ref_ev: baseline_e_ref_out,
1999 warnings: degenerate_normalization_warning(config)
2000 .into_iter()
2001 .collect(),
2002 })
2003}
2004
2005fn build_density_params(config: &UnifiedFitConfig) -> Vec<FitParameter> {
2008 config
2009 .initial_densities
2010 .iter()
2011 .enumerate()
2012 .map(|(i, &d)| {
2013 let name = config
2014 .isotope_names
2015 .get(i)
2016 .cloned()
2017 .unwrap_or_else(|| format!("isotope_{i}"));
2018 if config.density_is_fixed(i) {
2023 FitParameter::fixed(name, d)
2024 } else {
2025 FitParameter::non_negative(name, d)
2026 }
2027 })
2028 .collect()
2029}
2030
2031fn append_temperature_param(
2037 param_vec: &mut Vec<FitParameter>,
2038 config: &UnifiedFitConfig,
2039) -> Option<usize> {
2040 if !config.fit_temperature {
2041 return None;
2042 }
2043 let idx = param_vec.len();
2044 param_vec.push(FitParameter {
2045 name: "temperature_k".into(),
2046 value: config.temperature_k,
2047 lower: 1.0,
2048 upper: 5000.0,
2049 fixed: false,
2050 });
2051 Some(idx)
2052}
2053
2054fn detect_transmission_dips(measured: &[f64], energies: &[f64]) -> Vec<(f64, f64)> {
2060 let n = measured.len();
2061 if n < 5 || energies.len() != n {
2062 return Vec::new();
2063 }
2064 let mut sm = measured.to_vec();
2065 for i in 1..n - 1 {
2066 sm[i] = (measured[i - 1] + measured[i] + measured[i + 1]) / 3.0;
2067 }
2068 let mut sorted = sm.clone();
2070 sorted.sort_by(f64::total_cmp);
2071 let baseline = sorted[((0.9 * (n as f64 - 1.0)).round() as usize).min(n - 1)];
2072 let mut dips: Vec<(f64, f64)> = Vec::new();
2073 let mut max_depth = 0.0f64;
2074 for i in 1..n - 1 {
2075 if sm[i] < sm[i - 1] && sm[i] <= sm[i + 1] {
2076 let depth = baseline - sm[i];
2077 if depth > 0.0 {
2078 dips.push((energies[i], depth));
2079 max_depth = max_depth.max(depth);
2080 }
2081 }
2082 }
2083 dips.retain(|&(_, d)| d >= 0.2 * max_depth);
2086 dips
2087}
2088
2089fn peak_match_energy_scale_seed(
2105 measured: &[f64],
2106 energies: &[f64],
2107 config: &UnifiedFitConfig,
2108 flight_path_m: f64,
2109 t0_bounds: (f64, f64),
2110 l_scale_bounds: (f64, f64),
2111) -> Option<(f64, f64)> {
2112 let n = energies.len();
2113 if n < 5 || measured.len() != n {
2114 return None;
2115 }
2116 let refs: Vec<&ResonanceData> = config.resonance_data().iter().collect();
2118 let (e_lo, e_hi) = (
2119 energies[0].min(energies[n - 1]),
2120 energies[0].max(energies[n - 1]),
2121 );
2122 let res_e: Vec<f64> = nereids_physics::transmission::resonance_center_energies(&refs)
2123 .into_iter()
2124 .filter(|&e| e > e_lo && e < e_hi)
2125 .collect();
2126 if res_e.len() < 2 {
2127 return None;
2128 }
2129 let dips = detect_transmission_dips(measured, energies);
2130 if dips.len() < 2 {
2131 return None;
2132 }
2133 let min_spacing = res_e
2148 .windows(2)
2149 .map(|w| (w[1] - w[0]).abs())
2150 .fold(f64::INFINITY, f64::min);
2151 let mut grid_steps: Vec<f64> = energies.windows(2).map(|w| (w[1] - w[0]).abs()).collect();
2152 grid_steps.sort_by(f64::total_cmp);
2153 let grid_res = grid_steps.get(grid_steps.len() / 2).copied().unwrap_or(0.0);
2154 let match_tol = (0.5 * min_spacing).max(grid_res);
2155 let mut best: Vec<Option<(f64, f64)>> = vec![None; res_e.len()];
2156 for &(e_dip, depth) in &dips {
2157 let Some((k, &re)) = res_e
2158 .iter()
2159 .enumerate()
2160 .min_by(|(_, a), (_, b)| (*a - e_dip).abs().total_cmp(&(*b - e_dip).abs()))
2161 else {
2162 continue;
2163 };
2164 if (re - e_dip).abs() > match_tol {
2165 continue;
2166 }
2167 match best[k] {
2168 Some((_, d)) if d >= depth => {}
2169 _ => best[k] = Some((e_dip, depth)),
2170 }
2171 }
2172 let tof_factor = (0.5 * NEUTRON_MASS_KG / EV_TO_JOULES).sqrt() * 1.0e6;
2176 let c = tof_factor * flight_path_m;
2177 let pairs: Vec<(f64, f64)> = res_e
2178 .iter()
2179 .zip(best.iter())
2180 .filter_map(|(&re, b)| b.map(|(e_dip, _)| (c / re.sqrt(), c / e_dip.sqrt())))
2181 .collect();
2182 if pairs.len() < 2 {
2183 return None;
2184 }
2185 let m = pairs.len() as f64;
2187 let mean_x = pairs.iter().map(|p| p.0).sum::<f64>() / m;
2188 let mean_y = pairs.iter().map(|p| p.1).sum::<f64>() / m;
2189 let mut sxx = 0.0f64;
2190 let mut sxy = 0.0f64;
2191 for &(x, y) in &pairs {
2192 sxx += (x - mean_x) * (x - mean_x);
2193 sxy += (x - mean_x) * (y - mean_y);
2194 }
2195 if sxx <= 0.0 {
2196 return None;
2197 }
2198 let l_scale = sxy / sxx;
2199 let t0 = mean_y - l_scale * mean_x;
2200 if !t0.is_finite() || !l_scale.is_finite() {
2201 return None;
2202 }
2203 if t0 < t0_bounds.0
2208 || t0 > t0_bounds.1
2209 || l_scale < l_scale_bounds.0
2210 || l_scale > l_scale_bounds.1
2211 {
2212 return None;
2213 }
2214 Some((t0, l_scale))
2215}
2216
2217fn seed_energy_scale_in_params(
2221 param_vec: &mut [FitParameter],
2222 energy_scale_indices: Option<(usize, usize)>,
2223 measured_transmission: &[f64],
2224 config: &UnifiedFitConfig,
2225) {
2226 if !config.energy_scale_seed_enabled {
2229 return;
2230 }
2231 let Some((t0_idx, ls_idx)) = energy_scale_indices else {
2232 return;
2233 };
2234 let t0_b = (param_vec[t0_idx].lower, param_vec[t0_idx].upper);
2235 let ls_b = (param_vec[ls_idx].lower, param_vec[ls_idx].upper);
2236 if let Some((t0_seed, ls_seed)) = peak_match_energy_scale_seed(
2237 measured_transmission,
2238 config.energies(),
2239 config,
2240 config.flight_path_m,
2241 t0_b,
2242 ls_b,
2243 ) {
2244 param_vec[t0_idx].value = t0_seed;
2245 param_vec[ls_idx].value = ls_seed;
2246 }
2247}
2248
2249const ENERGY_SCALE_T0_BOUND_US: f64 = 10.0;
2253const ENERGY_SCALE_L_SCALE_LO: f64 = 0.99;
2257const ENERGY_SCALE_L_SCALE_HI: f64 = 1.01;
2259
2260fn append_energy_scale_params(
2266 param_vec: &mut Vec<FitParameter>,
2267 config: &UnifiedFitConfig,
2268) -> Option<(usize, usize)> {
2269 if !config.fit_energy_scale {
2270 return None;
2271 }
2272 let t0_idx = param_vec.len();
2273 param_vec.push(FitParameter {
2274 name: "t0_us".into(),
2275 value: config.t0_init_us,
2276 lower: -ENERGY_SCALE_T0_BOUND_US,
2277 upper: ENERGY_SCALE_T0_BOUND_US,
2278 fixed: false,
2279 });
2280 let ls_idx = param_vec.len();
2281 param_vec.push(FitParameter {
2282 name: "l_scale".into(),
2283 value: config.l_scale_init,
2284 lower: ENERGY_SCALE_L_SCALE_LO,
2285 upper: ENERGY_SCALE_L_SCALE_HI,
2286 fixed: false,
2287 });
2288 Some((t0_idx, ls_idx))
2289}
2290
2291pub(crate) fn validate_transmission_background(bg: &BackgroundConfig) -> Result<(), PipelineError> {
2306 if bg.fit_back_d != bg.fit_back_f {
2307 return Err(PipelineError::InvalidParameter(format!(
2308 "transmission_background: fit_back_d ({}) and fit_back_f ({}) \
2309 must both be true or both be false. The exponential tail \
2310 wrapper (BackD · exp(−BackF / √E)) requires both parameters \
2311 together; enabling only one leaves the other registered but \
2312 unused, silently producing the initial value as the fitted \
2313 result. Either enable both (to fit the exponential tail) or \
2314 disable both (4-term wrapper without exponential).",
2315 bg.fit_back_d, bg.fit_back_f,
2316 )));
2317 }
2318 Ok(())
2319}
2320
2321pub(crate) fn validate_multiplicative_baseline(
2337 config: &UnifiedFitConfig,
2338) -> Result<(), PipelineError> {
2339 let Some(bl) = config.multiplicative_baseline() else {
2340 return Ok(());
2341 };
2342 for (name, init, (lo, hi)) in [
2343 ("b0", bl.b0_init, bl.b0_bounds),
2344 ("b1", bl.b1_init, bl.b1_bounds),
2345 ("b2", bl.b2_init, bl.b2_bounds),
2346 ] {
2347 if !init.is_finite() {
2348 return Err(PipelineError::InvalidParameter(format!(
2349 "multiplicative baseline: {name}_init must be finite, got {init}"
2350 )));
2351 }
2352 if !lo.is_finite() || !hi.is_finite() || lo >= hi {
2353 return Err(PipelineError::InvalidParameter(format!(
2354 "multiplicative baseline: {name}_bounds must be finite with \
2355 lower < upper, got ({lo}, {hi})"
2356 )));
2357 }
2358 if init < lo || init > hi {
2359 return Err(PipelineError::InvalidParameter(format!(
2360 "multiplicative baseline: {name}_init = {init} lies outside \
2361 {name}_bounds ({lo}, {hi})"
2362 )));
2363 }
2364 }
2365 let e_ref = config.baseline_reference_energy();
2378 if !e_ref.is_finite() || e_ref <= 0.0 {
2379 return Err(PipelineError::InvalidParameter(format!(
2380 "multiplicative baseline: reference energy sqrt(E_min*E_max) is \
2381 invalid ({e_ref}) — check the energy grid"
2382 )));
2383 }
2384 let active_mask = nereids_fitting::active_mask::build_active_mask(
2385 config.energies(),
2386 config.fit_energy_range(),
2387 );
2388 for (i, &e) in config.energies().iter().enumerate() {
2389 if active_mask.as_ref().is_some_and(|m| !m[i]) {
2390 continue;
2391 }
2392 let z = (e / e_ref).ln();
2393 let b = bl.b0_init + bl.b1_init * z + bl.b2_init * z * z;
2394 let positive = b.is_finite() && b > 0.0;
2395 if !positive {
2396 return Err(PipelineError::InvalidParameter(format!(
2397 "multiplicative baseline: initial B(E) = {b} is not strictly \
2398 positive at E = {e} eV (inside the fit window) — adjust the \
2399 b0/b1/b2 inits"
2400 )));
2401 }
2402 }
2403 if config
2404 .transmission_background()
2405 .is_some_and(|bg| bg.fit_anorm)
2406 {
2407 return Err(PipelineError::InvalidParameter(
2408 "multiplicative baseline and a FREE SAMMY Anorm cannot be fitted \
2409 together: b0 and Anorm are degenerate normalizations. Set \
2410 BackgroundConfig::fit_anorm = false to combine the additive ABC \
2411 background with the baseline."
2412 .into(),
2413 ));
2414 }
2415 Ok(())
2416}
2417
2418pub(crate) fn degenerate_normalization_warning(config: &UnifiedFitConfig) -> Option<String> {
2424 let anorm_free = config
2425 .transmission_background()
2426 .is_some_and(|bg| bg.fit_anorm);
2427 if anorm_free && config.fit_temperature && config.n_free_density_params() >= 1 {
2428 Some(
2429 "fit configuration frees Anorm, temperature, AND at least one \
2430 density together — a degenerate normalization trio on real data \
2431 (observed: T ran to 4471 K with chi2/nu 932 and no warning). \
2432 Consider with_multiplicative_baseline (bounded normalization) \
2433 and/or with_fix_densities (known areal density)."
2434 .to_string(),
2435 )
2436 } else {
2437 None
2438 }
2439}
2440
2441pub(crate) fn validate_precomputed_cross_sections(
2468 config: &UnifiedFitConfig,
2469) -> Result<(), PipelineError> {
2470 let Some(xs) = config.precomputed_cross_sections() else {
2471 return Ok(());
2472 };
2473 if xs.is_empty() {
2474 return Err(PipelineError::ShapeMismatch(
2475 "precomputed_cross_sections must not be empty".into(),
2476 ));
2477 }
2478
2479 let n_e = config.energies().len();
2480 for (i, row) in xs.iter().enumerate() {
2481 if row.len() != n_e {
2482 return Err(PipelineError::ShapeMismatch(format!(
2483 "precomputed_cross_sections row {i} has length {} but config.energies \
2484 has {n_e}",
2485 row.len(),
2486 )));
2487 }
2488 if let Some(j) = row.iter().position(|s| !s.is_finite()) {
2495 return Err(PipelineError::ShapeMismatch(format!(
2496 "precomputed_cross_sections row {i} has non-finite σ at energy index {j}: {}",
2497 row[j],
2498 )));
2499 }
2500 }
2501
2502 let n_params = config.n_density_params();
2503 let member_rows = config.density_indices.as_ref().map(|di| di.len());
2506 let row_ok = xs.len() == n_params || member_rows == Some(xs.len());
2507 if !row_ok {
2508 let expected = match member_rows {
2509 Some(m) if m != n_params => format!("{n_params} (collapsed) or {m} (per-member)"),
2510 _ => format!("{n_params}"),
2511 };
2512 return Err(PipelineError::ShapeMismatch(format!(
2513 "precomputed_cross_sections has {} rows but expected {expected}",
2514 xs.len(),
2515 )));
2516 }
2517
2518 if let Some((work_xs, layout)) = &config.precomputed_work_cross_sections {
2531 let n_work = layout.energies.len();
2532 if work_xs.is_empty() {
2533 return Err(PipelineError::ShapeMismatch(
2534 "precomputed_work_cross_sections must not be empty".into(),
2535 ));
2536 }
2537 for (i, row) in work_xs.iter().enumerate() {
2538 if row.len() != n_work {
2539 return Err(PipelineError::ShapeMismatch(format!(
2540 "precomputed_work_cross_sections row {i} has length {} but \
2541 work_layout has {n_work} working-grid energies",
2542 row.len(),
2543 )));
2544 }
2545 if let Some(j) = row.iter().position(|s| !s.is_finite()) {
2546 return Err(PipelineError::ShapeMismatch(format!(
2547 "precomputed_work_cross_sections row {i} has non-finite σ at \
2548 working-grid index {j}: {}",
2549 row[j],
2550 )));
2551 }
2552 }
2553 if work_xs.len() != xs.len() {
2556 return Err(PipelineError::ShapeMismatch(format!(
2557 "precomputed_work_cross_sections has {} rows but \
2558 precomputed_cross_sections has {} — both index the same density \
2559 mapping and must agree",
2560 work_xs.len(),
2561 xs.len(),
2562 )));
2563 }
2564 if layout.data_indices.len() != n_e {
2568 return Err(PipelineError::ShapeMismatch(format!(
2569 "precomputed_work_cross_sections layout maps {} data points but \
2570 config.energies has {n_e}",
2571 layout.data_indices.len(),
2572 )));
2573 }
2574 if let Some(&bad) = layout.data_indices.iter().find(|&&idx| idx >= n_work) {
2575 return Err(PipelineError::ShapeMismatch(format!(
2576 "precomputed_work_cross_sections layout index {bad} is out of \
2577 range for {n_work} working-grid energies",
2578 )));
2579 }
2580 }
2581 Ok(())
2582}
2583
2584pub(crate) fn count_free_params(config: &UnifiedFitConfig) -> usize {
2614 let mut n_free = config.n_free_density_params();
2618 if config.fit_temperature {
2619 n_free += 1;
2620 }
2621 if config.fit_energy_scale {
2622 n_free += 2;
2623 }
2624 if let Some(bg) = config.transmission_background.as_ref() {
2625 n_free += usize::from(bg.fit_anorm);
2626 n_free += usize::from(bg.fit_back_a);
2627 n_free += usize::from(bg.fit_back_b);
2628 n_free += usize::from(bg.fit_back_c);
2629 n_free += usize::from(bg.fit_back_d);
2630 n_free += usize::from(bg.fit_back_f);
2631 }
2632 if let Some(bl) = config.multiplicative_baseline.as_ref() {
2633 n_free += usize::from(bl.fit_b0);
2634 n_free += usize::from(bl.fit_b1);
2635 n_free += usize::from(bl.fit_b2);
2636 }
2637 n_free
2638}
2639
2640pub(crate) fn required_active_bins(config: &UnifiedFitConfig) -> usize {
2659 count_free_params(config).max(2)
2660}
2661
2662fn append_background_params(
2663 param_vec: &mut Vec<FitParameter>,
2664 bg: &BackgroundConfig,
2665) -> BackgroundIndices {
2666 let anorm = param_vec.len();
2671 param_vec.push(if bg.fit_anorm {
2672 FitParameter {
2673 name: "anorm".into(),
2674 value: bg.anorm_init,
2675 lower: 0.5,
2676 upper: 2.0,
2677 fixed: false,
2678 }
2679 } else {
2680 FitParameter::fixed("anorm", bg.anorm_init)
2681 });
2682 let back_a = param_vec.len();
2688 param_vec.push(if bg.fit_back_a {
2689 FitParameter {
2690 name: "back_a".into(),
2691 value: bg.back_a_init,
2692 lower: -0.5,
2693 upper: 0.5,
2694 fixed: false,
2695 }
2696 } else {
2697 FitParameter::fixed("back_a", bg.back_a_init)
2698 });
2699 let back_b = param_vec.len();
2700 param_vec.push(if bg.fit_back_b {
2701 FitParameter {
2702 name: "back_b".into(),
2703 value: bg.back_b_init,
2704 lower: -0.5,
2705 upper: 0.5,
2706 fixed: false,
2707 }
2708 } else {
2709 FitParameter::fixed("back_b", bg.back_b_init)
2710 });
2711 let back_c = param_vec.len();
2712 param_vec.push(if bg.fit_back_c {
2713 FitParameter {
2714 name: "back_c".into(),
2715 value: bg.back_c_init,
2716 lower: -0.5,
2717 upper: 0.5,
2718 fixed: false,
2719 }
2720 } else {
2721 FitParameter::fixed("back_c", bg.back_c_init)
2722 });
2723
2724 let back_d = if bg.fit_back_d {
2732 let idx = param_vec.len();
2733 param_vec.push(FitParameter {
2734 name: "back_d".into(),
2735 value: bg.back_d_init,
2736 lower: 0.0,
2737 upper: 1.0,
2738 fixed: false,
2739 });
2740 Some(idx)
2741 } else {
2742 None
2743 };
2744 let back_f = if bg.fit_back_f {
2745 let idx = param_vec.len();
2746 param_vec.push(FitParameter {
2747 name: "back_f".into(),
2748 value: bg.back_f_init,
2749 lower: 0.0,
2750 upper: 100.0,
2751 fixed: false,
2752 });
2753 Some(idx)
2754 } else {
2755 None
2756 };
2757
2758 BackgroundIndices {
2759 anorm,
2760 back_a,
2761 back_b,
2762 back_c,
2763 back_d,
2764 back_f,
2765 }
2766}
2767
2768fn append_multiplicative_baseline_params(
2780 param_vec: &mut Vec<FitParameter>,
2781 bl: &MultiplicativeBaselineConfig,
2782) -> BaselineIndices {
2783 let b0 = param_vec.len();
2784 param_vec.push(if bl.fit_b0 {
2785 FitParameter {
2786 name: "baseline_b0".into(),
2787 value: bl.b0_init,
2788 lower: bl.b0_bounds.0,
2789 upper: bl.b0_bounds.1,
2790 fixed: false,
2791 }
2792 } else {
2793 FitParameter::fixed("baseline_b0", bl.b0_init)
2794 });
2795 let b1 = param_vec.len();
2796 param_vec.push(if bl.fit_b1 {
2797 FitParameter {
2798 name: "baseline_b1".into(),
2799 value: bl.b1_init,
2800 lower: bl.b1_bounds.0,
2801 upper: bl.b1_bounds.1,
2802 fixed: false,
2803 }
2804 } else {
2805 FitParameter::fixed("baseline_b1", bl.b1_init)
2806 });
2807 let b2 = param_vec.len();
2808 param_vec.push(if bl.fit_b2 {
2809 FitParameter {
2810 name: "baseline_b2".into(),
2811 value: bl.b2_init,
2812 lower: bl.b2_bounds.0,
2813 upper: bl.b2_bounds.1,
2814 fixed: false,
2815 }
2816 } else {
2817 FitParameter::fixed("baseline_b2", bl.b2_init)
2818 });
2819
2820 BaselineIndices { b0, b1, b2 }
2821}
2822
2823fn build_energy_scale_transmission_model(
2847 config: &UnifiedFitConfig,
2848 t0_idx: usize,
2849 ls_idx: usize,
2850 temperature_index: Option<usize>,
2851) -> Result<Box<dyn FitModel>, PipelineError> {
2852 let instrument = config
2853 .resolution
2854 .clone()
2855 .map(|r| Arc::new(InstrumentParams { resolution: r }));
2856 let n_iso = config.resonance_data.len();
2862 let density_indices = config
2863 .density_indices
2864 .clone()
2865 .unwrap_or_else(|| (0..n_iso).collect());
2866 let density_ratios = config
2867 .density_ratios
2868 .clone()
2869 .unwrap_or_else(|| vec![1.0; n_iso]);
2870 let mut es_model = EnergyScaleTransmissionModel::new(
2871 Arc::new(config.resonance_data.clone()),
2872 Arc::new(density_indices),
2873 Arc::new(density_ratios),
2874 config.temperature_k,
2875 config.energies.clone(),
2876 config.flight_path_m,
2877 t0_idx,
2878 ls_idx,
2879 instrument,
2880 )
2881 .with_temperature_index(temperature_index)
2885 .map_err(|e| PipelineError::InvalidParameter(format!("energy-scale model: {e}")))?;
2886 if let Some(method) = config.tzero_jacobian_method {
2887 es_model = es_model.with_jacobian_method(method);
2888 }
2889 Ok(Box::new(es_model))
2890}
2891
2892fn build_transmission_model(
2894 config: &UnifiedFitConfig,
2895 n_density_params: usize,
2896 temperature_index: Option<usize>,
2897) -> Result<Box<dyn FitModel>, PipelineError> {
2898 let n_params = config.n_density_params();
2899
2900 let computed_xs_storage;
2917 let effective_precomputed: Option<&Arc<Vec<Vec<f64>>>> = if config.fit_temperature {
2918 None
2919 } else if let Some(xs) = &config.precomputed_cross_sections {
2920 Some(xs)
2921 } else {
2922 let instrument = config
2923 .resolution
2924 .clone()
2925 .map(|r| Arc::new(InstrumentParams { resolution: r }));
2926 let working = nereids_physics::transmission::broadened_cross_sections_on_working_grid(
2927 config.energies(),
2928 &config.resonance_data,
2929 config.temperature_k,
2930 instrument.as_deref(),
2931 None,
2932 )
2933 .map_err(PipelineError::Transmission)?;
2934 if working.layout.is_identity() {
2935 computed_xs_storage = Arc::new(working.sigma);
2937 Some(&computed_xs_storage)
2938 } else {
2939 let (density_indices, density_ratios) = (
2945 config
2946 .density_indices
2947 .clone()
2948 .unwrap_or_else(|| (0..n_params).collect()),
2949 config
2950 .density_ratios
2951 .clone()
2952 .unwrap_or_else(|| vec![1.0; config.resonance_data.len()]),
2953 );
2954 if density_indices.len() != working.sigma.len()
2960 || density_ratios.len() != working.sigma.len()
2961 {
2962 return Err(PipelineError::InvalidParameter(format!(
2963 "density mapping length mismatch: {} density_indices / {} \
2964 density_ratios for {} per-isotope cross-section rows",
2965 density_indices.len(),
2966 density_ratios.len(),
2967 working.sigma.len(),
2968 )));
2969 }
2970 let n_e = working.sigma[0].len();
2971 let mut eff = vec![vec![0.0f64; n_e]; n_params];
2972 for ((&idx, &ratio), member_xs) in density_indices
2973 .iter()
2974 .zip(density_ratios.iter())
2975 .zip(working.sigma.iter())
2976 {
2977 for (j, &sigma) in member_xs.iter().enumerate() {
2978 eff[idx][j] += ratio * sigma;
2979 }
2980 }
2981 return Ok(Box::new(PrecomputedTransmissionModel {
2982 cross_sections: Arc::new(eff),
2983 density_indices: Arc::new((0..n_params).collect()),
2984 energies: instrument
2985 .as_ref()
2986 .map(|_| Arc::new(config.energies.clone())),
2987 instrument,
2988 resolution_plan: None,
2989 sparse_cubature_plan: None,
2990 sparse_scalar_plan: None,
2991 work_layout: Some(Arc::new(working.layout)),
2992 }));
2993 }
2994 };
2995
2996 if !config.fit_temperature
2997 && let Some(xs) = effective_precomputed
2998 {
2999 let collapse_by_groups = |xs: &Arc<Vec<Vec<f64>>>| -> Arc<Vec<Vec<f64>>> {
3004 if let (Some(di), Some(dr)) = (&config.density_indices, &config.density_ratios)
3005 && xs.len() == di.len()
3006 && di.len() == dr.len()
3007 {
3008 let n_e = xs[0].len();
3009 let mut eff = vec![vec![0.0f64; n_e]; n_params];
3010 for ((&idx, &ratio), member_xs) in di.iter().zip(dr.iter()).zip(xs.iter()) {
3011 for (j, &sigma) in member_xs.iter().enumerate() {
3012 eff[idx][j] += ratio * sigma;
3013 }
3014 }
3015 return Arc::new(eff);
3016 }
3017 Arc::clone(xs)
3018 };
3019
3020 let (effective_xs, work_layout): (
3028 Arc<Vec<Vec<f64>>>,
3029 Option<Arc<nereids_physics::transmission::WorkingGridLayout>>,
3030 ) = match &config.precomputed_work_cross_sections {
3031 Some((work_xs, layout)) => (collapse_by_groups(work_xs), Some(Arc::clone(layout))),
3032 None => (collapse_by_groups(xs), None),
3033 };
3034 let instrument = config
3037 .resolution
3038 .clone()
3039 .map(|r| Arc::new(InstrumentParams { resolution: r }));
3040 let resolution_plan = if instrument.is_some() {
3041 config.precomputed_resolution_plan.clone()
3042 } else {
3043 None
3044 };
3045 let sparse_cubature_plan = if instrument.is_some() {
3046 config.precomputed_sparse_cubature_plan.clone()
3047 } else {
3048 None
3049 };
3050 let sparse_scalar_plan = if instrument.is_some() {
3051 config.precomputed_sparse_scalar_plan.clone()
3052 } else {
3053 None
3054 };
3055 return Ok(Box::new(PrecomputedTransmissionModel {
3056 cross_sections: effective_xs,
3057 density_indices: Arc::new((0..n_params).collect()),
3058 energies: instrument
3059 .as_ref()
3060 .map(|_| Arc::new(config.energies.clone())),
3061 instrument,
3062 resolution_plan,
3063 sparse_cubature_plan,
3064 sparse_scalar_plan,
3065 work_layout,
3066 }));
3067 }
3068
3069 let instrument = config
3070 .resolution
3071 .clone()
3072 .map(|r| Arc::new(InstrumentParams { resolution: r }));
3073
3074 let base_xs = config.precomputed_base_xs.clone();
3075 let density_ratios = config
3076 .density_ratios
3077 .clone()
3078 .unwrap_or_else(|| vec![1.0; n_density_params]);
3079 let density_indices = config
3080 .density_indices
3081 .clone()
3082 .unwrap_or_else(|| (0..n_density_params).collect());
3083 let resolution_plan = if instrument.is_some() {
3084 config.precomputed_resolution_plan.clone()
3085 } else {
3086 None
3087 };
3088 let sparse_cubature_plan = if instrument.is_some() {
3089 config.precomputed_sparse_cubature_plan.clone()
3090 } else {
3091 None
3092 };
3093 let sparse_scalar_plan = if instrument.is_some() {
3094 config.precomputed_sparse_scalar_plan.clone()
3095 } else {
3096 None
3097 };
3098 Ok(Box::new(
3099 TransmissionFitModel::new(
3100 config.energies.clone(),
3101 config.resonance_data.clone(),
3102 config.temperature_k,
3103 instrument,
3104 (density_indices, density_ratios),
3105 temperature_index,
3106 base_xs,
3107 )?
3108 .with_resolution_plan(resolution_plan)
3109 .with_sparse_cubature_plan(sparse_cubature_plan)
3110 .with_sparse_scalar_plan(sparse_scalar_plan),
3111 ))
3112}
3113
3114fn poisson_to_lm_result(
3116 model: &dyn FitModel,
3117 measured_t: &[f64],
3118 sigma: &[f64],
3119 pr: &poisson::PoissonResult,
3120 params: &ParameterSet,
3121 scale_by_chi2: bool,
3122) -> Result<LmResult, PipelineError> {
3123 let n_free = params.n_free();
3124 let dof = measured_t.len().saturating_sub(n_free).max(1);
3125 let y_model = model.evaluate(&pr.params)?;
3126 let chi_sq: f64 = measured_t
3127 .iter()
3128 .zip(y_model.iter())
3129 .zip(sigma.iter())
3130 .map(|((obs, mdl), s)| {
3131 let residual = obs - mdl;
3132 (residual * residual) / (s * s).max(1e-30)
3133 })
3134 .sum();
3135 let reduced_chi_squared = chi_sq / dof as f64;
3136
3137 let (covariance, uncertainties) =
3152 if scale_by_chi2 && reduced_chi_squared.is_finite() && reduced_chi_squared > 0.0 {
3153 let sigma_factor = reduced_chi_squared.sqrt();
3154 let covariance = pr.covariance.as_ref().map(|cov| {
3155 let mut scaled = cov.clone();
3156 for v in scaled.data.iter_mut() {
3157 *v *= reduced_chi_squared;
3158 }
3159 scaled
3160 });
3161 let uncertainties = pr
3162 .uncertainties
3163 .as_ref()
3164 .map(|unc| unc.iter().map(|s| s * sigma_factor).collect());
3165 (covariance, uncertainties)
3166 } else {
3167 (pr.covariance.clone(), pr.uncertainties.clone())
3168 };
3169
3170 Ok(LmResult {
3171 chi_squared: chi_sq,
3172 reduced_chi_squared,
3173 iterations: pr.iterations,
3174 converged: pr.converged,
3175 params: pr.params.clone(),
3176 covariance,
3177 uncertainties,
3178 })
3179}
3180
3181fn free_uncertainty(free_indices: &[usize], unc_all: &[f64], full_index: usize) -> Option<f64> {
3193 free_indices
3194 .iter()
3195 .position(|&fi| fi == full_index)
3196 .and_then(|pos| unc_all.get(pos).copied())
3197}
3198
3199fn extract_result(
3201 config: &UnifiedFitConfig,
3202 result: &LmResult,
3203 n_density_params: usize,
3204 free_indices: &[usize],
3205 bg_indices: Option<BackgroundIndices>,
3206 bl_indices: Option<BaselineIndices>,
3207) -> Result<SpectrumFitResult, PipelineError> {
3208 let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();
3209
3210 let (baseline, baseline_e_ref_ev) = if let Some(bli) = bl_indices {
3214 (
3215 Some([
3216 result.params[bli.b0],
3217 result.params[bli.b1],
3218 result.params[bli.b2],
3219 ]),
3220 Some(config.baseline_reference_energy()),
3221 )
3222 } else {
3223 (None, None)
3224 };
3225
3226 let (anorm, background, back_d, back_f): (f64, [f64; 3], Option<f64>, Option<f64>) =
3227 if let Some(bi) = bg_indices {
3228 let bd = bi.back_d.map(|i| result.params[i]);
3231 let bf = bi.back_f.map(|i| result.params[i]);
3232 (
3233 result.params[bi.anorm],
3234 [
3235 result.params[bi.back_a],
3236 result.params[bi.back_b],
3237 result.params[bi.back_c],
3238 ],
3239 bd,
3240 bf,
3241 )
3242 } else {
3243 (1.0, [0.0, 0.0, 0.0], None, None)
3244 };
3245
3246 let (uncertainties, temperature_k, temperature_k_unc) = if result.converged {
3247 match &result.uncertainties {
3248 Some(unc_all) => {
3249 let (temp_k, temp_unc) = if config.fit_temperature {
3258 (
3259 Some(result.params[n_density_params]),
3260 Some(
3261 free_uncertainty(free_indices, unc_all, n_density_params)
3262 .unwrap_or(f64::NAN),
3263 ),
3264 )
3265 } else {
3266 (None, None)
3267 };
3268 let unc: Vec<f64> = (0..n_density_params)
3269 .map(|i| free_uncertainty(free_indices, unc_all, i).unwrap_or(f64::NAN))
3270 .collect();
3271 (Some(unc), temp_k, temp_unc)
3272 }
3273 None => {
3274 let temp_k = if config.fit_temperature {
3275 Some(result.params[n_density_params])
3276 } else {
3277 None
3278 };
3279 (None, temp_k, None)
3280 }
3281 }
3282 } else {
3283 let temp_k = if config.fit_temperature {
3284 Some(result.params[n_density_params])
3285 } else {
3286 None
3287 };
3288 (None, temp_k, None)
3289 };
3290
3291 Ok(SpectrumFitResult {
3292 densities,
3293 uncertainties,
3294 reduced_chi_squared: result.reduced_chi_squared,
3295 converged: result.converged,
3296 iterations: result.iterations,
3297 temperature_k,
3298 temperature_k_unc,
3299 anorm,
3300 background,
3301 back_d,
3302 back_f,
3303 t0_us: None,
3304 l_scale: None,
3305 energy_scale_flight_path_m: None,
3306 deviance_per_dof: None,
3307 baseline,
3308 baseline_e_ref_ev,
3309 warnings: degenerate_normalization_warning(config)
3310 .into_iter()
3311 .collect(),
3312 })
3313}
3314
3315pub struct ModelJacobianResult {
3323 pub jacobian: lm::FlatMatrix,
3325 pub fisher: lm::FlatMatrix,
3327 pub model_prediction: Vec<f64>,
3329 pub param_names: Vec<String>,
3331}
3332
3333pub fn evaluate_jacobian_and_fisher(
3360 config: &UnifiedFitConfig,
3361 flux: &[f64],
3362 background: &[f64],
3363) -> Result<ModelJacobianResult, PipelineError> {
3364 validate_precomputed_cross_sections(config)?;
3373
3374 if config.multiplicative_baseline.is_some() {
3382 return Err(PipelineError::InvalidParameter(
3383 "evaluate_jacobian_and_fisher does not support a multiplicative \
3384 baseline: its parameter layout is frozen for research-script \
3385 compatibility (kl_b0/kl_b1 background stand-ins). Remove \
3386 with_multiplicative_baseline from the config for this helper."
3387 .into(),
3388 ));
3389 }
3390
3391 let n_density_params = config.n_density_params();
3392
3393 let mut param_vec = build_density_params(config);
3400
3401 let temperature_index = if config.fit_temperature {
3402 let idx = param_vec.len();
3403 param_vec.push(FitParameter {
3404 name: "temperature_k".into(),
3405 value: config.temperature_k,
3406 lower: 1.0,
3407 upper: 5000.0,
3408 fixed: false,
3409 });
3410 Some(idx)
3411 } else {
3412 None
3413 };
3414
3415 let kl_bg = if config.transmission_background.is_some() {
3416 let base = param_vec.len();
3417 param_vec.push(FitParameter {
3418 name: "kl_b0".into(),
3419 value: 0.0,
3420 lower: 0.0,
3421 upper: 0.5,
3422 fixed: false,
3423 });
3424 param_vec.push(FitParameter {
3425 name: "kl_b1".into(),
3426 value: 0.0,
3427 lower: 0.0,
3428 upper: 0.5,
3429 fixed: false,
3430 });
3431 Some((base, base + 1))
3432 } else {
3433 None
3434 };
3435
3436 let counts_bg = if let Some(bg) = config.counts_background() {
3437 let alpha1_idx = param_vec.len();
3438 param_vec.push(if bg.fit_alpha_1 {
3439 FitParameter {
3440 name: "alpha_1".into(),
3441 value: bg.alpha_1_init,
3442 lower: 0.0,
3443 upper: 10.0,
3444 fixed: false,
3445 }
3446 } else {
3447 FitParameter::fixed("alpha_1", bg.alpha_1_init)
3448 });
3449 let alpha2_idx = param_vec.len();
3450 param_vec.push(if bg.fit_alpha_2 {
3451 FitParameter {
3452 name: "alpha_2".into(),
3453 value: bg.alpha_2_init,
3454 lower: 0.0,
3455 upper: 10.0,
3456 fixed: false,
3457 }
3458 } else {
3459 FitParameter::fixed("alpha_2", bg.alpha_2_init)
3460 });
3461 Some((alpha1_idx, alpha2_idx))
3462 } else {
3463 None
3464 };
3465
3466 let params = ParameterSet::new(param_vec);
3467 let all_vals = params.all_values();
3468 let free_idx = params.free_indices();
3469 let n_free = free_idx.len();
3470
3471 let param_names: Vec<String> = free_idx
3473 .iter()
3474 .map(|&i| params.params[i].name.to_string())
3475 .collect();
3476
3477 let config_with_xs;
3502 let effective_config =
3503 if config.fit_temperature || config.precomputed_work_cross_sections.is_some() {
3504 config
3505 } else {
3506 let instrument = config
3507 .resolution
3508 .clone()
3509 .map(|r| Arc::new(InstrumentParams { resolution: r }));
3510 let working = nereids_physics::transmission::broadened_cross_sections_on_working_grid(
3511 config.energies(),
3512 &config.resonance_data,
3513 config.temperature_k,
3514 instrument.as_deref(),
3515 None,
3516 )
3517 .map_err(PipelineError::Transmission)?;
3518 config_with_xs = if working.layout.is_identity() {
3519 config
3521 .clone()
3522 .with_precomputed_cross_sections(Arc::new(working.sigma))
3523 } else {
3524 let data_xs: Vec<Vec<f64>> = working
3529 .sigma
3530 .iter()
3531 .map(|s| working.layout.extract(s))
3532 .collect();
3533 config
3534 .clone()
3535 .with_precomputed_cross_sections(Arc::new(data_xs))
3536 .with_precomputed_work_cross_sections(
3537 Arc::new(working.sigma),
3538 Arc::new(working.layout),
3539 )
3540 };
3541 &config_with_xs
3542 };
3543
3544 let t_model = build_transmission_model(effective_config, n_density_params, temperature_index)?;
3546
3547 let evaluate_and_jacobian =
3550 |model: &dyn FitModel| -> Result<(Vec<f64>, lm::FlatMatrix), PipelineError> {
3551 let y_model = model.evaluate(&all_vals)?;
3552 let jac = model
3553 .analytical_jacobian(&all_vals, &free_idx, &y_model)
3554 .ok_or_else(|| {
3555 PipelineError::InvalidParameter(
3556 "analytical Jacobian not available for this model configuration".into(),
3557 )
3558 })?;
3559 Ok((y_model, jac))
3560 };
3561
3562 let (y_model, jac) = if let Some((b0_idx, b1_idx)) = kl_bg {
3563 let inv_sqrt_e: Vec<f64> = config
3564 .energies()
3565 .iter()
3566 .map(|&e| 1.0 / e.max(1e-10).sqrt())
3567 .collect();
3568 let wrapped = poisson::TransmissionKLBackgroundModel {
3569 inner: &*t_model,
3570 inv_sqrt_energies: inv_sqrt_e,
3571 b0_index: b0_idx,
3572 b1_index: b1_idx,
3573 n_params: params.params.len(),
3574 };
3575 if let Some((a1, a2)) = counts_bg {
3576 let cm = poisson::CountsBackgroundScaleModel {
3577 transmission_model: &wrapped,
3578 flux,
3579 background,
3580 alpha1_index: a1,
3581 alpha2_index: a2,
3582 n_params: params.params.len(),
3583 };
3584 evaluate_and_jacobian(&cm)?
3585 } else {
3586 let cm = poisson::CountsModel {
3587 transmission_model: &wrapped,
3588 flux,
3589 background,
3590 n_params: params.params.len(),
3591 };
3592 evaluate_and_jacobian(&cm)?
3593 }
3594 } else if let Some((a1, a2)) = counts_bg {
3595 let cm = poisson::CountsBackgroundScaleModel {
3596 transmission_model: &*t_model,
3597 flux,
3598 background,
3599 alpha1_index: a1,
3600 alpha2_index: a2,
3601 n_params: params.params.len(),
3602 };
3603 evaluate_and_jacobian(&cm)?
3604 } else {
3605 let cm = poisson::CountsModel {
3606 transmission_model: &*t_model,
3607 flux,
3608 background,
3609 n_params: params.params.len(),
3610 };
3611 evaluate_and_jacobian(&cm)?
3612 };
3613
3614 let mut fisher = lm::FlatMatrix::zeros(n_free, n_free);
3616 for (i, &mu_i) in y_model.iter().enumerate() {
3617 let mu_inv = 1.0 / mu_i.max(1e-30);
3618 for a in 0..n_free {
3619 let ja = jac.get(i, a);
3620 for b in 0..=a {
3621 let jb = jac.get(i, b);
3622 *fisher.get_mut(a, b) += ja * jb * mu_inv;
3623 if a != b {
3624 *fisher.get_mut(b, a) += ja * jb * mu_inv;
3625 }
3626 }
3627 }
3628 }
3629
3630 Ok(ModelJacobianResult {
3631 jacobian: jac,
3632 fisher,
3633 model_prediction: y_model,
3634 param_names,
3635 })
3636}
3637
3638#[derive(Debug, PartialEq)]
3642pub enum FitConfigError {
3643 EmptyEnergies,
3645 EmptyResonanceData,
3647 DensityCountMismatch { densities: usize, isotopes: usize },
3649 NameCountMismatch { names: usize, isotopes: usize },
3651 GroupMemberCountMismatch {
3653 group_name: String,
3654 rd_count: usize,
3655 member_count: usize,
3656 },
3657 GroupMemberIsotopeMismatch {
3659 group_name: String,
3660 expected_z: u32,
3661 expected_a: u32,
3662 got_z: u32,
3663 got_a: u32,
3664 },
3665 NonFiniteTemperature(f64),
3667 NegativeTemperature(f64),
3669 InvalidFitEnergyRange(&'static str),
3671 DensityFreezeBeforeGroups,
3676}
3677
3678impl fmt::Display for FitConfigError {
3679 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
3680 match self {
3681 Self::EmptyEnergies => write!(f, "energy grid must be non-empty"),
3682 Self::EmptyResonanceData => write!(f, "resonance_data must be non-empty"),
3683 Self::DensityCountMismatch {
3684 densities,
3685 isotopes,
3686 } => write!(
3687 f,
3688 "initial_densities length ({densities}) must match number of density parameters ({isotopes})"
3689 ),
3690 Self::NameCountMismatch { names, isotopes } => write!(
3691 f,
3692 "isotope_names length ({names}) must match resonance_data length ({isotopes})"
3693 ),
3694 Self::GroupMemberCountMismatch {
3695 group_name,
3696 rd_count,
3697 member_count,
3698 } => write!(
3699 f,
3700 "group '{group_name}': provided {rd_count} ResonanceData but group has {member_count} members"
3701 ),
3702 Self::GroupMemberIsotopeMismatch {
3703 group_name,
3704 expected_z,
3705 expected_a,
3706 got_z,
3707 got_a,
3708 } => write!(
3709 f,
3710 "group '{group_name}': expected Z={expected_z} A={expected_a} but got Z={got_z} A={got_a}"
3711 ),
3712 Self::NonFiniteTemperature(v) => {
3713 write!(f, "temperature must be finite, got {v}")
3714 }
3715 Self::NegativeTemperature(v) => {
3716 write!(f, "temperature must be non-negative, got {v}")
3717 }
3718 Self::InvalidFitEnergyRange(msg) => {
3719 write!(f, "invalid fit_energy_range: {msg}")
3720 }
3721 Self::DensityFreezeBeforeGroups => write!(
3722 f,
3723 "density freeze (with_fix_densities / with_density_free) must be configured \
3724 after with_groups: grouping redefines the density parameters"
3725 ),
3726 }
3727 }
3728}
3729
3730impl std::error::Error for FitConfigError {}
3731
3732#[derive(Debug, Clone)]
3734pub struct SpectrumFitResult {
3735 pub densities: Vec<f64>,
3737 pub uncertainties: Option<Vec<f64>>,
3741 pub reduced_chi_squared: f64,
3743 pub converged: bool,
3745 pub iterations: usize,
3747 pub temperature_k: Option<f64>,
3749 pub temperature_k_unc: Option<f64>,
3762 pub anorm: f64,
3766 pub background: [f64; 3],
3779 pub back_d: Option<f64>,
3787 pub back_f: Option<f64>,
3794 pub t0_us: Option<f64>,
3797 pub l_scale: Option<f64>,
3800 pub energy_scale_flight_path_m: Option<f64>,
3807 pub deviance_per_dof: Option<f64>,
3816 pub baseline: Option<[f64; 3]>,
3828 pub baseline_e_ref_ev: Option<f64>,
3835 pub warnings: Vec<String>,
3844}
3845
3846impl SpectrumFitResult {
3847 pub fn corrected_energies(
3870 &self,
3871 nominal_energies: &[f64],
3872 ) -> Option<Result<Vec<f64>, PipelineError>> {
3873 match (self.t0_us, self.l_scale, self.energy_scale_flight_path_m) {
3874 (Some(t0), Some(l_scale), Some(flight_path_m)) => Some(
3875 nereids_fitting::resolution_calib::corrected_energy_grid(
3876 nominal_energies,
3877 t0,
3878 l_scale,
3879 flight_path_m,
3880 )
3881 .map_err(|e| PipelineError::InvalidParameter(e.to_string())),
3882 ),
3883 _ => None,
3884 }
3885 }
3886}
3887
3888#[cfg(test)]
3889mod tests {
3890 use super::*;
3891 use nereids_endf::resonance::test_support::{
3892 hf178_mlbw_two_resonances, synthetic_single_resonance, u238_single_resonance,
3893 u238_three_resonances,
3894 };
3895 use nereids_fitting::lm::FitModel;
3896 use nereids_fitting::transmission_model::{
3897 EnergyScaleJacobianMethod, EnergyScaleTransmissionModel,
3898 };
3899 use nereids_physics::transmission as phys_transmission;
3900
3901 fn invert_dense(a: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
3905 let n = a.len();
3906 let mut m: Vec<Vec<f64>> = a.to_vec();
3907 let mut inv: Vec<Vec<f64>> = (0..n)
3908 .map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect())
3909 .collect();
3910 for col in 0..n {
3911 let mut piv = col;
3913 for r in (col + 1)..n {
3914 if m[r][col].abs() > m[piv][col].abs() {
3915 piv = r;
3916 }
3917 }
3918 if m[piv][col].abs() < 1e-300 {
3919 return None;
3920 }
3921 m.swap(col, piv);
3922 inv.swap(col, piv);
3923 let d = m[col][col];
3924 for j in 0..n {
3925 m[col][j] /= d;
3926 inv[col][j] /= d;
3927 }
3928 for r in 0..n {
3929 if r == col {
3930 continue;
3931 }
3932 let f = m[r][col];
3933 for j in 0..n {
3934 m[r][j] -= f * m[col][j];
3935 inv[r][j] -= f * inv[col][j];
3936 }
3937 }
3938 }
3939 Some(inv)
3940 }
3941
3942 #[test]
3945 fn detect_transmission_dips_finds_clear_dips() {
3946 let energies: Vec<f64> = (0..120).map(|i| 1.0 + (i as f64) * 0.5).collect();
3947 let mut t = vec![1.0_f64; 120];
3948 for v in t.iter_mut().take(33).skip(28) {
3950 *v = 0.4;
3951 }
3952 for v in t.iter_mut().take(83).skip(78) {
3953 *v = 0.6;
3954 }
3955 let dips = detect_transmission_dips(&t, &energies);
3956 assert_eq!(dips.len(), 2, "expected 2 dips, got {dips:?}");
3957 let mut de: Vec<f64> = dips.iter().map(|&(e, _)| e).collect();
3958 de.sort_by(f64::total_cmp);
3959 assert!(
3961 (de[0] - energies[30]).abs() <= 1.0,
3962 "first dip at {} eV",
3963 de[0]
3964 );
3965 assert!(
3966 (de[1] - energies[80]).abs() <= 1.0,
3967 "second dip at {} eV",
3968 de[1]
3969 );
3970 }
3971
3972 #[test]
3978 fn baseline_reference_energy_honours_fit_energy_range() {
3979 let data = hf178_mlbw_two_resonances();
3980 let mut energies: Vec<f64> = (0..400).map(|i| 8.0 + (i as f64) * 0.0925).collect();
3982 energies.push(3211.0 * 3211.0 / 8.0); let config = UnifiedFitConfig::new(
3984 energies.clone(),
3985 vec![data],
3986 vec!["Hf-178".into()],
3987 293.6,
3988 None,
3989 vec![0.1],
3990 )
3991 .unwrap();
3992 let full = nereids_fitting::transmission_model::baseline_reference_energy(&energies);
3993 assert!((config.baseline_reference_energy() - full).abs() < 1e-6);
3995 let windowed = config.with_fit_energy_range(Some((8.0, 45.0))).unwrap();
3997 let e_ref = windowed.baseline_reference_energy();
3998 assert!(
3999 e_ref > 8.0 && e_ref < 45.0,
4000 "windowed E_ref = {e_ref} eV must lie inside the 8–45 eV fit window"
4001 );
4002 assert!(
4003 (e_ref - full).abs() > 100.0,
4004 "windowed E_ref must differ from the buggy full-grid value {full}"
4005 );
4006 }
4007
4008 #[test]
4013 fn peak_match_energy_scale_seed_identity() {
4014 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
4016 let (t_obs, _sigma) = synthetic_transmission(&data, 0.1, &energies);
4017 let config = UnifiedFitConfig::new(
4018 energies.clone(),
4019 vec![data],
4020 vec!["Hf-178".into()],
4021 293.6,
4022 None,
4023 vec![0.1],
4024 )
4025 .unwrap()
4026 .with_energy_scale(0.0, 1.0, 25.0);
4027 let (t0, l_scale) = peak_match_energy_scale_seed(
4028 &t_obs,
4029 config.energies(),
4030 &config,
4031 25.0,
4032 (-10.0, 10.0),
4033 (0.99, 1.01),
4034 )
4035 .expect("seed should be Some with 2 resonances and detectable dips");
4036 assert!(
4038 t0.abs() < 0.5,
4039 "t0 should be ≈0 for un-shifted data, got {t0}"
4040 );
4041 assert!(
4042 (l_scale - 1.0).abs() < 5e-3,
4043 "L_scale should be ≈1 for un-shifted data, got {l_scale}"
4044 );
4045 }
4046
4047 #[test]
4052 fn peak_match_energy_scale_seed_recovers_nonidentity() {
4053 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
4057 let density = 0.1_f64;
4058 let flight_path = 25.0_f64;
4059 let (t0_true, l_scale_true) = (2.0_f64, 1.006_f64);
4060 let model = EnergyScaleTransmissionModel::new(
4063 Arc::new(vec![data.clone()]),
4064 Arc::new(vec![0]),
4065 Arc::new(vec![1.0]),
4066 293.6,
4067 energies.clone(),
4068 flight_path,
4069 1, 2, None,
4072 );
4073 let t_obs = model.evaluate(&[density, t0_true, l_scale_true]).unwrap();
4074 let config = UnifiedFitConfig::new(
4075 energies.clone(),
4076 vec![data],
4077 vec!["Hf-178".into()],
4078 293.6,
4079 None,
4080 vec![density],
4081 )
4082 .unwrap()
4083 .with_energy_scale(0.0, 1.0, flight_path);
4084 let (t0, l_scale) = peak_match_energy_scale_seed(
4085 &t_obs,
4086 config.energies(),
4087 &config,
4088 flight_path,
4089 (-10.0, 10.0),
4090 (0.99, 1.01),
4091 )
4092 .expect("seed should be Some for a clean shifted two-resonance spectrum");
4093 assert!(
4097 (t0 - t0_true).abs() < 0.6,
4098 "seed should recover t0 ≈ {t0_true}, got {t0}"
4099 );
4100 assert!(
4101 (l_scale - l_scale_true).abs() < 4e-3,
4102 "seed should recover L_scale ≈ {l_scale_true}, got {l_scale}"
4103 );
4104 assert!(
4106 (t0 - t0_true).abs() < (0.0 - t0_true).abs()
4107 && (l_scale - l_scale_true).abs() < (1.0 - l_scale_true).abs(),
4108 "seed must improve on the cold start"
4109 );
4110 }
4111
4112 #[test]
4117 fn peak_match_energy_scale_seed_none_on_featureless_spectrum() {
4118 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
4120 let t_obs: Vec<f64> = energies.iter().map(|&e| (-5.0 / e.sqrt()).exp()).collect();
4122 let config = UnifiedFitConfig::new(
4123 energies.clone(),
4124 vec![data],
4125 vec!["Hf-178".into()],
4126 293.6,
4127 None,
4128 vec![0.1],
4129 )
4130 .unwrap()
4131 .with_energy_scale(0.0, 1.0, 25.0);
4132 assert!(
4133 peak_match_energy_scale_seed(
4134 &t_obs,
4135 config.energies(),
4136 &config,
4137 25.0,
4138 (-10.0, 10.0),
4139 (0.99, 1.01),
4140 )
4141 .is_none(),
4142 "a featureless heavily-absorbing spectrum has no resonance dips ⇒ \
4143 seed must return None (cold-start fallback)"
4144 );
4145 }
4146
4147 #[test]
4154 fn peak_match_energy_scale_seed_handles_duplicate_resonance_energies() {
4155 use nereids_physics::transmission::{SampleParams, forward_model};
4156
4157 let iso_a = synthetic_single_resonance(72, 178, 176.0, 7.8);
4159 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();
4162 let density = 0.05_f64;
4163 let sample = SampleParams::new(
4164 293.6,
4165 vec![
4166 (iso_a.clone(), density),
4167 (iso_b.clone(), density),
4168 (iso_c.clone(), density),
4169 ],
4170 )
4171 .unwrap();
4172 let t_obs = forward_model(&energies, &sample, None).unwrap();
4173 let config = UnifiedFitConfig::new(
4174 energies.clone(),
4175 vec![iso_a, iso_b, iso_c],
4176 vec!["A".into(), "B".into(), "C".into()],
4177 293.6,
4178 None,
4179 vec![density, density, density],
4180 )
4181 .unwrap()
4182 .with_energy_scale(0.0, 1.0, 25.0);
4183 let (t0, l_scale) = peak_match_energy_scale_seed(
4184 &t_obs,
4185 config.energies(),
4186 &config,
4187 25.0,
4188 (-10.0, 10.0),
4189 (0.99, 1.01),
4190 )
4191 .expect(
4192 "duplicate resonance energies must not collapse match_tol to 0 — \
4193 seed should be Some (was silently None pre-fix)",
4194 );
4195 assert!(t0.abs() < 0.5, "identity calibration: t0 ≈ 0, got {t0}");
4196 assert!(
4197 (l_scale - 1.0).abs() < 5e-3,
4198 "identity calibration: L_scale ≈ 1, got {l_scale}"
4199 );
4200 }
4201
4202 #[test]
4205 fn test_input_data_transmission_n_energies() {
4206 let data = InputData::Transmission {
4207 transmission: vec![0.9, 0.8, 0.7],
4208 uncertainty: vec![0.01, 0.01, 0.01],
4209 };
4210 assert_eq!(data.n_energies(), 3);
4211 assert!(!data.is_counts());
4212 }
4213
4214 #[test]
4215 fn test_input_data_counts_n_energies() {
4216 let data = InputData::Counts {
4217 sample_counts: vec![10.0, 20.0, 30.0, 40.0],
4218 open_beam_counts: vec![100.0, 100.0, 100.0, 100.0],
4219 };
4220 assert_eq!(data.n_energies(), 4);
4221 assert!(data.is_counts());
4222 }
4223
4224 #[test]
4225 fn test_input_data_counts_with_nuisance() {
4226 let data = InputData::CountsWithNuisance {
4227 sample_counts: vec![5.0, 6.0],
4228 flux: vec![100.0, 100.0],
4229 background: vec![0.5, 0.5],
4230 };
4231 assert_eq!(data.n_energies(), 2);
4232 assert!(data.is_counts());
4233 }
4234
4235 #[test]
4236 fn test_solver_config_default_is_auto() {
4237 let cfg = SolverConfig::default();
4238 assert!(matches!(cfg, SolverConfig::Auto));
4239 }
4240
4241 #[test]
4242 fn test_counts_background_config_default() {
4243 let cfg = CountsBackgroundConfig::default();
4244 assert!((cfg.alpha_1_init - 1.0).abs() < f64::EPSILON);
4245 assert!((cfg.alpha_2_init - 1.0).abs() < f64::EPSILON);
4246 assert!(!cfg.fit_alpha_1);
4247 assert!(!cfg.fit_alpha_2);
4248 }
4249
4250 fn synthetic_transmission(
4254 data: &ResonanceData,
4255 true_density: f64,
4256 energies: &[f64],
4257 ) -> (Vec<f64>, Vec<f64>) {
4258 let model = PrecomputedTransmissionModel {
4259 cross_sections: Arc::new(vec![
4260 phys_transmission::broadened_cross_sections(
4261 energies,
4262 std::slice::from_ref(data),
4263 0.0,
4264 None,
4265 None,
4266 )
4267 .unwrap()
4268 .into_iter()
4269 .next()
4270 .unwrap(),
4271 ]),
4272 density_indices: Arc::new(vec![0]),
4273 energies: None,
4274 instrument: None,
4275 resolution_plan: None,
4276 sparse_cubature_plan: None,
4277 sparse_scalar_plan: None,
4278 work_layout: None,
4279 };
4280 let t = model.evaluate(&[true_density]).unwrap();
4281 let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4282 (t, sigma)
4283 }
4284
4285 fn synthetic_counts(
4287 data: &ResonanceData,
4288 true_density: f64,
4289 energies: &[f64],
4290 i0: f64,
4291 ) -> (Vec<f64>, Vec<f64>) {
4292 let (t, _) = synthetic_transmission(data, true_density, energies);
4293 let open_beam: Vec<f64> = vec![i0; energies.len()];
4294 let sample: Vec<f64> = t.iter().map(|&v| (v * i0).round().max(0.0)).collect();
4295 (sample, open_beam)
4296 }
4297
4298 #[test]
4299 fn test_typed_transmission_lm_recovers_density() {
4300 let data = u238_single_resonance();
4301 let true_density = 0.002;
4302 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4303 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
4304
4305 let config = UnifiedFitConfig::new(
4306 energies,
4307 vec![data],
4308 vec!["U-238".into()],
4309 0.0,
4310 None,
4311 vec![0.001],
4312 )
4313 .unwrap()
4314 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4315
4316 let input = InputData::Transmission {
4317 transmission: t,
4318 uncertainty: sigma,
4319 };
4320
4321 let result = fit_spectrum_typed(&input, &config).unwrap();
4322 assert!(result.converged, "LM should converge");
4323 let fitted = result.densities[0];
4324 assert!(
4325 (fitted - true_density).abs() / true_density < 0.05,
4326 "density: fitted={fitted}, true={true_density}"
4327 );
4328 }
4329
4330 fn precomputed_xs_fixture() -> (UnifiedFitConfig, InputData) {
4333 let data = u238_single_resonance();
4334 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
4335 let (t, sigma) = synthetic_transmission(&data, 0.001, &energies);
4336 let config = UnifiedFitConfig::new(
4337 energies,
4338 vec![data],
4339 vec!["U-238".into()],
4340 0.0,
4341 None,
4342 vec![0.001],
4343 )
4344 .unwrap()
4345 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4346 let input = InputData::Transmission {
4347 transmission: t,
4348 uncertainty: sigma,
4349 };
4350 (config, input)
4351 }
4352
4353 #[test]
4354 fn test_fit_rejects_empty_precomputed_cross_sections() {
4355 let (config, input) = precomputed_xs_fixture();
4358 let config = config.with_precomputed_cross_sections(Arc::new(Vec::new()));
4359 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4360 assert!(
4361 matches!(err, PipelineError::ShapeMismatch(_)),
4362 "expected ShapeMismatch, got {err:?}"
4363 );
4364 assert!(err.to_string().contains("must not be empty"));
4365 }
4366
4367 #[test]
4368 fn test_fit_rejects_wrong_row_count_precomputed_cross_sections() {
4369 let (config, input) = precomputed_xs_fixture();
4371 let n_e = config.energies().len();
4372 let bad = vec![vec![1.0; n_e], vec![1.0; n_e]];
4373 let config = config.with_precomputed_cross_sections(Arc::new(bad));
4374 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4375 assert!(
4376 matches!(err, PipelineError::ShapeMismatch(_)),
4377 "expected ShapeMismatch, got {err:?}"
4378 );
4379 assert!(err.to_string().contains("rows"));
4380 }
4381
4382 #[test]
4383 fn test_fit_rejects_wrong_energy_length_precomputed_cross_sections() {
4384 let (config, input) = precomputed_xs_fixture();
4386 let n_e = config.energies().len();
4387 let bad = vec![vec![1.0; n_e + 3]];
4388 let config = config.with_precomputed_cross_sections(Arc::new(bad));
4389 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4390 assert!(
4391 matches!(err, PipelineError::ShapeMismatch(_)),
4392 "expected ShapeMismatch, got {err:?}"
4393 );
4394 assert!(err.to_string().contains("config.energies"));
4395 }
4396
4397 #[test]
4398 fn test_fit_accepts_correct_precomputed_cross_sections() {
4399 let (config, input) = precomputed_xs_fixture();
4402 let n_e = config.energies().len();
4403 let xs = phys_transmission::broadened_cross_sections(
4404 config.energies(),
4405 config.resonance_data(),
4406 0.0,
4407 None,
4408 None,
4409 )
4410 .unwrap();
4411 assert_eq!(xs.len(), 1);
4412 assert_eq!(xs[0].len(), n_e);
4413 let config = config.with_precomputed_cross_sections(Arc::new(xs));
4414 let result = fit_spectrum_typed(&input, &config);
4417 assert!(
4418 result.is_ok(),
4419 "correctly-shaped precomputed XS must pass validation, got {result:?}"
4420 );
4421 }
4422
4423 #[test]
4424 fn test_fit_rejects_non_finite_precomputed_cross_sections() {
4425 let (config, input) = precomputed_xs_fixture();
4431 let n_e = config.energies().len();
4432 let bad = vec![vec![f64::NAN; n_e]];
4433 let config = config.with_precomputed_cross_sections(Arc::new(bad));
4434 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4435 assert!(
4436 matches!(err, PipelineError::ShapeMismatch(_)),
4437 "expected ShapeMismatch, got {err:?}"
4438 );
4439 assert!(
4440 err.to_string().contains("non-finite"),
4441 "error should mention non-finite σ, got: {err}"
4442 );
4443
4444 let (config, input) = precomputed_xs_fixture();
4446 let mut row = vec![1.0; n_e];
4447 row[n_e / 2] = f64::INFINITY;
4448 let config = config.with_precomputed_cross_sections(Arc::new(vec![row]));
4449 let err = fit_spectrum_typed(&input, &config).unwrap_err();
4450 assert!(
4451 matches!(err, PipelineError::ShapeMismatch(_)),
4452 "expected ShapeMismatch for +inf σ, got {err:?}"
4453 );
4454 }
4455
4456 #[test]
4457 fn test_evaluate_jacobian_rejects_malformed_precomputed_cross_sections() {
4458 let (config, _input) = precomputed_xs_fixture();
4463 let n_e = config.energies().len();
4464 let flux = vec![1.0; n_e];
4465 let background = vec![0.0; n_e];
4466
4467 let empty = config
4470 .clone()
4471 .with_precomputed_cross_sections(Arc::new(Vec::new()));
4472 match evaluate_jacobian_and_fisher(&empty, &flux, &background) {
4473 Err(PipelineError::ShapeMismatch(msg)) => {
4474 assert!(msg.contains("must not be empty"), "got: {msg}");
4475 }
4476 Err(other) => panic!("expected ShapeMismatch for empty XS, got {other:?}"),
4477 Ok(_) => panic!("empty precomputed XS must be rejected"),
4478 }
4479
4480 let nan = config.with_precomputed_cross_sections(Arc::new(vec![vec![f64::NAN; n_e]]));
4482 match evaluate_jacobian_and_fisher(&nan, &flux, &background) {
4483 Err(PipelineError::ShapeMismatch(msg)) => {
4484 assert!(msg.contains("non-finite"), "got: {msg}");
4485 }
4486 Err(other) => panic!("expected ShapeMismatch for NaN σ, got {other:?}"),
4487 Ok(_) => panic!("non-finite precomputed σ must be rejected"),
4488 }
4489 }
4490
4491 #[test]
4492 fn test_extract_result_drops_uncertainties_when_unconverged() {
4493 let data = u238_single_resonance();
4494 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
4495 let config = UnifiedFitConfig::new(
4496 energies,
4497 vec![data],
4498 vec!["U-238".into()],
4499 293.6,
4500 None,
4501 vec![0.001],
4502 )
4503 .unwrap();
4504
4505 let result = LmResult {
4506 chi_squared: 1.0,
4507 reduced_chi_squared: 1.0,
4508 iterations: 5,
4509 converged: false,
4510 params: vec![0.001],
4511 covariance: Some(lm::FlatMatrix::zeros(1, 1)),
4512 uncertainties: Some(vec![0.123]),
4513 };
4514
4515 let extracted = extract_result(&config, &result, 1, &[0], None, None).unwrap();
4517 assert!(!extracted.converged);
4518 assert!(
4519 extracted.uncertainties.is_none(),
4520 "pipeline must not surface uncertainties from an unconverged fit"
4521 );
4522 assert!(extracted.temperature_k_unc.is_none());
4523 }
4524
4525 #[test]
4526 fn test_typed_counts_kl_recovers_density() {
4527 let data = u238_single_resonance();
4528 let true_density = 0.002;
4529 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4530 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4531
4532 let config = UnifiedFitConfig::new(
4533 energies,
4534 vec![data],
4535 vec!["U-238".into()],
4536 0.0,
4537 None,
4538 vec![0.001],
4539 )
4540 .unwrap()
4541 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4542
4543 let input = InputData::Counts {
4544 sample_counts: sample,
4545 open_beam_counts: open_beam,
4546 };
4547
4548 let result = fit_spectrum_typed(&input, &config).unwrap();
4549 assert!(result.converged, "KL on counts should converge");
4550 let fitted = result.densities[0];
4551 assert!(
4552 (fitted - true_density).abs() / true_density < 0.10,
4553 "density: fitted={fitted}, true={true_density}"
4554 );
4555 }
4556
4557 #[test]
4558 fn test_typed_counts_kl_low_counts_recovers_density() {
4559 let data = u238_single_resonance();
4561 let true_density = 0.0005;
4562 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4563 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 10.0);
4564
4565 let config = UnifiedFitConfig::new(
4566 energies,
4567 vec![data],
4568 vec!["U-238".into()],
4569 0.0,
4570 None,
4571 vec![0.001],
4572 )
4573 .unwrap()
4574 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4575
4576 let input = InputData::Counts {
4577 sample_counts: sample,
4578 open_beam_counts: open_beam,
4579 };
4580
4581 let result = fit_spectrum_typed(&input, &config).unwrap();
4582 assert!(result.converged, "KL on low counts should converge");
4583 let fitted = result.densities[0];
4584 assert!(
4586 (fitted - true_density).abs() / true_density < 0.30,
4587 "density: fitted={fitted}, true={true_density}"
4588 );
4589 }
4590
4591 #[test]
4592 fn test_typed_transmission_kl_recovers_density() {
4593 let data = u238_single_resonance();
4594 let true_density = 0.0005;
4595 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4596 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
4597
4598 let config = UnifiedFitConfig::new(
4599 energies,
4600 vec![data],
4601 vec!["U-238".into()],
4602 0.0,
4603 None,
4604 vec![0.001],
4605 )
4606 .unwrap()
4607 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
4608
4609 let input = InputData::Transmission {
4610 transmission: t,
4611 uncertainty: sigma,
4612 };
4613
4614 let result = fit_spectrum_typed(&input, &config).unwrap();
4615 assert!(result.converged, "KL on transmission should converge");
4616 let fitted = result.densities[0];
4617 assert!(
4618 (fitted - true_density).abs() / true_density < 0.05,
4619 "density: fitted={fitted}, true={true_density}"
4620 );
4621 }
4622
4623 #[test]
4624 fn test_typed_counts_lm_auto_converts() {
4625 let data = u238_single_resonance();
4627 let true_density = 0.0005;
4628 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4629 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4630
4631 let config = UnifiedFitConfig::new(
4632 energies,
4633 vec![data],
4634 vec!["U-238".into()],
4635 0.0,
4636 None,
4637 vec![0.001],
4638 )
4639 .unwrap()
4640 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4641
4642 let input = InputData::Counts {
4643 sample_counts: sample,
4644 open_beam_counts: open_beam,
4645 };
4646
4647 let result = fit_spectrum_typed(&input, &config).unwrap();
4648 assert!(
4649 result.converged,
4650 "LM on auto-converted counts should converge"
4651 );
4652 let fitted = result.densities[0];
4653 assert!(
4654 (fitted - true_density).abs() / true_density < 0.10,
4655 "density: fitted={fitted}, true={true_density}"
4656 );
4657 }
4658
4659 #[test]
4660 fn test_typed_auto_solver_selects_kl_for_counts() {
4661 let data = u238_single_resonance();
4662 let true_density = 0.0005;
4663 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4664 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
4665
4666 let config = UnifiedFitConfig::new(
4668 energies,
4669 vec![data],
4670 vec!["U-238".into()],
4671 0.0,
4672 None,
4673 vec![0.001],
4674 )
4675 .unwrap(); let input = InputData::Counts {
4678 sample_counts: sample,
4679 open_beam_counts: open_beam,
4680 };
4681
4682 let result = fit_spectrum_typed(&input, &config).unwrap();
4683 assert!(
4684 result.converged,
4685 "Auto solver on counts should use KL and converge"
4686 );
4687 }
4688
4689 #[test]
4690 fn test_typed_transmission_with_background_lm() {
4691 let data = u238_single_resonance();
4692 let true_density = 0.0005;
4693 let true_anorm = 0.95;
4694 let true_back_a = 0.02;
4695 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4696
4697 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
4699 let t_bg: Vec<f64> = t_pure
4700 .iter()
4701 .map(|&v| true_anorm * v + true_back_a)
4702 .collect();
4703 let sigma: Vec<f64> = t_bg.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4704
4705 let config = UnifiedFitConfig::new(
4706 energies,
4707 vec![data],
4708 vec!["U-238".into()],
4709 0.0,
4710 None,
4711 vec![0.001],
4712 )
4713 .unwrap()
4714 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
4715 max_iter: 500,
4716 ..LmConfig::default()
4717 }))
4718 .with_transmission_background(BackgroundConfig::default());
4719
4720 let input = InputData::Transmission {
4721 transmission: t_bg,
4722 uncertainty: sigma,
4723 };
4724
4725 let result = fit_spectrum_typed(&input, &config).unwrap();
4726 assert!(
4727 result.converged,
4728 "LM+BG should converge on noiseless synthetic data (chi2r={}, iter={})",
4729 result.reduced_chi_squared, result.iterations
4730 );
4731 assert!(
4732 (result.densities[0] - true_density).abs() / true_density < 0.05,
4733 "density: fitted={}, true={true_density}",
4734 result.densities[0]
4735 );
4736 assert!(
4737 (result.anorm - true_anorm).abs() / true_anorm < 0.05,
4738 "anorm: fitted={}, true={true_anorm}",
4739 result.anorm
4740 );
4741 }
4742
4743 #[test]
4744 fn test_typed_counts_with_nuisance_rejects_lm() {
4745 let data = u238_single_resonance();
4746 let energies: Vec<f64> = (0..10).map(|i| 1.0 + (i as f64) * 0.5).collect();
4747
4748 let config = UnifiedFitConfig::new(
4749 energies,
4750 vec![data],
4751 vec!["U-238".into()],
4752 0.0,
4753 None,
4754 vec![0.001],
4755 )
4756 .unwrap()
4757 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4758
4759 let input = InputData::CountsWithNuisance {
4760 sample_counts: vec![10.0; 10],
4761 flux: vec![100.0; 10],
4762 background: vec![0.0; 10],
4763 };
4764
4765 let result = fit_spectrum_typed(&input, &config);
4766 assert!(result.is_err(), "CountsWithNuisance + LM should error");
4767 }
4768
4769 fn seeded_gaussian(n: usize, seed: u64) -> Vec<f64> {
4773 let mut state = seed | 1;
4774 let mut unif = || {
4775 state = state
4776 .wrapping_mul(6364136223846793005)
4777 .wrapping_add(1442695040888963407);
4778 (((state >> 11) as f64) / ((1u64 << 53) as f64)).clamp(1e-12, 1.0 - 1e-12)
4779 };
4780 (0..n)
4781 .map(|_| {
4782 let u1 = unif();
4783 let u2 = unif();
4784 (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos()
4785 })
4786 .collect()
4787 }
4788
4789 fn synthetic_transmission_at_temp(
4790 data: &ResonanceData,
4791 true_density: f64,
4792 temperature_k: f64,
4793 energies: &[f64],
4794 ) -> (Vec<f64>, Vec<f64>) {
4795 let xs = phys_transmission::broadened_cross_sections(
4796 energies,
4797 std::slice::from_ref(data),
4798 temperature_k,
4799 None,
4800 None,
4801 )
4802 .unwrap();
4803 let model = PrecomputedTransmissionModel {
4804 cross_sections: Arc::new(xs),
4805 density_indices: Arc::new(vec![0]),
4806 energies: None,
4807 instrument: None,
4808 resolution_plan: None,
4809 sparse_cubature_plan: None,
4810 sparse_scalar_plan: None,
4811 work_layout: None,
4812 };
4813 let t = model.evaluate(&[true_density]).unwrap();
4814 let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
4815 (t, sigma)
4816 }
4817
4818 #[test]
4819 fn test_typed_poisson_kl_with_temperature() {
4820 let data = u238_single_resonance();
4821 let true_density = 0.0005;
4822 let true_temp = 350.0;
4823 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4824 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
4825
4826 let config = UnifiedFitConfig::new(
4827 energies,
4828 vec![data],
4829 vec!["U-238".into()],
4830 300.0, None,
4832 vec![0.001],
4833 )
4834 .unwrap()
4835 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
4836 max_iter: 500,
4837 ..PoissonConfig::default()
4838 }))
4839 .with_fit_temperature(true);
4840
4841 let input = InputData::Transmission {
4842 transmission: t,
4843 uncertainty: sigma,
4844 };
4845
4846 let result = fit_spectrum_typed(&input, &config).unwrap();
4847
4848 let fitted_density = result.densities[0];
4850 assert!(
4851 (fitted_density - true_density).abs() / true_density < 0.01,
4852 "density: fitted={fitted_density}, true={true_density}, ratio={}",
4853 (fitted_density - true_density).abs() / true_density,
4854 );
4855
4856 let fitted_temp = result
4858 .temperature_k
4859 .expect("temperature_k should be Some when fit_temperature=true");
4860 assert!(
4861 (fitted_temp - true_temp).abs() < 1.0,
4862 "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
4863 (fitted_temp - true_temp).abs(),
4864 );
4865 }
4866
4867 #[test]
4868 fn test_typed_poisson_kl_with_temperature_and_background() {
4869 let data = u238_single_resonance();
4870 let true_density = 0.0005;
4871 let true_temp = 350.0;
4872 let true_b0 = 0.012;
4873 let true_b1 = 0.008;
4874 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
4875 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
4876 let measured_t: Vec<f64> = t
4877 .iter()
4878 .zip(energies.iter())
4879 .map(|(&ti, &e)| ti + true_b0 + true_b1 / e.sqrt())
4880 .collect();
4881
4882 let config = UnifiedFitConfig::new(
4883 energies,
4884 vec![data],
4885 vec!["U-238".into()],
4886 300.0,
4887 None,
4888 vec![0.001],
4889 )
4890 .unwrap()
4891 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
4892 max_iter: 120,
4893 gauss_newton_lambda: 1e-4,
4894 ..PoissonConfig::default()
4895 }))
4896 .with_fit_temperature(true)
4897 .with_transmission_background(BackgroundConfig::default());
4898
4899 let input = InputData::Transmission {
4900 transmission: measured_t,
4901 uncertainty: sigma,
4902 };
4903
4904 let result = fit_spectrum_typed(&input, &config).unwrap();
4905
4906 assert!(result.converged, "fit did not converge: {result:?}");
4907 assert!(
4908 result.iterations <= 80,
4909 "expected KL background+temperature fit to converge well before max_iter; got {}",
4910 result.iterations,
4911 );
4912
4913 let fitted_density = result.densities[0];
4914 assert!(
4915 (fitted_density - true_density).abs() / true_density < 0.02,
4916 "density: fitted={fitted_density}, true={true_density}, ratio={}",
4917 (fitted_density - true_density).abs() / true_density,
4918 );
4919
4920 let fitted_temp = result
4921 .temperature_k
4922 .expect("temperature_k should be Some when fit_temperature=true");
4923 assert!(
4924 (fitted_temp - true_temp).abs() < 3.0,
4925 "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
4926 (fitted_temp - true_temp).abs(),
4927 );
4928
4929 assert!(
4930 (result.background[0] - true_b0).abs() < 5e-3,
4931 "background b0: fitted={}, true={}",
4932 result.background[0],
4933 true_b0,
4934 );
4935 assert!(
4936 (result.background[1] - true_b1).abs() < 5e-3,
4937 "background b1: fitted={}, true={}",
4938 result.background[1],
4939 true_b1,
4940 );
4941 }
4942
4943 #[test]
4947 fn test_grouped_fit_spectrum_round_trip() {
4948 use nereids_core::types::IsotopeGroup;
4949
4950 let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
4952 let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);
4953
4954 let iso1 = nereids_core::types::Isotope::new(92, 235).unwrap();
4956 let iso2 = nereids_core::types::Isotope::new(92, 238).unwrap();
4957 let group =
4958 IsotopeGroup::custom("U (60/40)".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();
4959
4960 let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
4961 let true_density = 0.0005;
4962
4963 let sample = nereids_physics::transmission::SampleParams::new(
4965 0.0,
4966 vec![
4967 (rd1.clone(), true_density * 0.6),
4968 (rd2.clone(), true_density * 0.4),
4969 ],
4970 )
4971 .unwrap();
4972 let transmission =
4973 nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
4974 let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();
4975
4976 let config = UnifiedFitConfig::new(
4978 energies.clone(),
4979 vec![rd1.clone()],
4980 vec!["placeholder".into()],
4981 0.0,
4982 None,
4983 vec![0.001],
4984 )
4985 .unwrap()
4986 .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
4987 .unwrap()
4988 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
4989
4990 let input = InputData::Transmission {
4991 transmission,
4992 uncertainty,
4993 };
4994
4995 let result = fit_spectrum_typed(&input, &config).unwrap();
4996
4997 assert_eq!(result.densities.len(), 1, "should have 1 group density");
4999 let fitted = result.densities[0];
5000 let rel_error = (fitted - true_density).abs() / true_density;
5001 assert!(
5002 rel_error < 0.01,
5003 "group density: fitted={fitted}, true={true_density}, rel_error={rel_error}"
5004 );
5005 assert!(result.converged, "fit should converge");
5006 }
5007
5008 #[test]
5009 fn test_grouped_poisson_kl_with_temperature_and_background_noiseless() {
5010 use nereids_core::types::IsotopeGroup;
5011
5012 let rd1 = synthetic_single_resonance(72, 176, 8.5, 5.0);
5013 let rd2 = synthetic_single_resonance(72, 178, 17.0, 7.5);
5014 let rd3 = synthetic_single_resonance(72, 180, 29.0, 6.0);
5015
5016 let hf176 = nereids_core::types::Isotope::new(72, 176).unwrap();
5017 let hf178 = nereids_core::types::Isotope::new(72, 178).unwrap();
5018 let hf180 = nereids_core::types::Isotope::new(72, 180).unwrap();
5019 let group = IsotopeGroup::custom(
5020 "Hf-like (3 member)".into(),
5021 vec![(hf176, 0.2), (hf178, 0.5), (hf180, 0.3)],
5022 )
5023 .unwrap();
5024
5025 let energies: Vec<f64> = (0..300).map(|i| 1.0 + (49.0 * i as f64) / 299.0).collect();
5026 let true_density = 0.001;
5027 let true_temp = 400.0;
5028 let true_b0 = 0.012;
5029 let true_b1 = 0.008;
5030
5031 let sample = nereids_physics::transmission::SampleParams::new(
5032 true_temp,
5033 vec![
5034 (rd1.clone(), true_density * 0.2),
5035 (rd2.clone(), true_density * 0.5),
5036 (rd3.clone(), true_density * 0.3),
5037 ],
5038 )
5039 .unwrap();
5040 let pure_t =
5041 nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
5042 let measured_t: Vec<f64> = pure_t
5043 .iter()
5044 .zip(energies.iter())
5045 .map(|(&t, &e)| t + true_b0 + true_b1 / e.sqrt())
5046 .collect();
5047 let sigma = vec![0.001; energies.len()];
5048
5049 let config = UnifiedFitConfig::new(
5050 energies.clone(),
5051 vec![rd1.clone()],
5052 vec!["placeholder".into()],
5053 293.6,
5054 None,
5055 vec![0.0008],
5056 )
5057 .unwrap()
5058 .with_groups(&[(&group, &[rd1, rd2, rd3])], vec![0.0008])
5059 .unwrap()
5060 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
5061 max_iter: 200,
5062 gauss_newton_lambda: 1e-4,
5063 ..PoissonConfig::default()
5064 }))
5065 .with_fit_temperature(true)
5066 .with_transmission_background(BackgroundConfig::default());
5067
5068 let input = InputData::Transmission {
5069 transmission: measured_t,
5070 uncertainty: sigma,
5071 };
5072
5073 let result = fit_spectrum_typed(&input, &config).unwrap();
5074
5075 assert!(result.converged, "fit did not converge: {result:?}");
5076 assert!(
5079 (result.densities[0] - true_density).abs() / true_density < 0.01,
5080 "density: fitted={}, true={true_density}",
5081 result.densities[0]
5082 );
5083 let fitted_temp = result
5084 .temperature_k
5085 .expect("temperature_k should be Some when fit_temperature=true");
5086 assert!(
5087 (fitted_temp - true_temp).abs() < 8.0,
5088 "temperature: fitted={fitted_temp}, true={true_temp}",
5089 );
5090 let e_mid: f64 = 10.0;
5094 let bg_total = (result.anorm - 1.0)
5095 + result.background[0]
5096 + result.background[1] / e_mid.sqrt()
5097 + result.background[2] * e_mid.sqrt();
5098 let true_bg_mid = true_b0 + true_b1 / e_mid.sqrt();
5099 assert!(
5100 (bg_total - true_bg_mid).abs() < 0.02,
5101 "total bg at E={e_mid}: fitted={bg_total:.6}, true={true_bg_mid:.6}",
5102 );
5103 }
5104
5105 #[test]
5108 fn test_kl_counts_returns_density_uncertainty() {
5109 let data = u238_single_resonance();
5110 let true_density = 0.002;
5111 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5112 let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);
5113
5114 let config = UnifiedFitConfig::new(
5115 energies,
5116 vec![data],
5117 vec!["U-238".into()],
5118 0.0,
5119 None,
5120 vec![0.001],
5121 )
5122 .unwrap()
5123 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
5124
5125 let input = InputData::Counts {
5126 sample_counts: sample,
5127 open_beam_counts: open_beam,
5128 };
5129 let result = fit_spectrum_typed(&input, &config).unwrap();
5130 assert!(result.converged);
5131 let unc = result
5132 .uncertainties
5133 .as_ref()
5134 .expect("KL 1D fit should return density uncertainties");
5135 assert_eq!(unc.len(), 1);
5136 assert!(
5137 unc[0].is_finite() && unc[0] > 0.0,
5138 "density unc = {}",
5139 unc[0]
5140 );
5141 assert!(
5142 unc[0] < result.densities[0],
5143 "unc ({}) should be < density ({}) for high-count data",
5144 unc[0],
5145 result.densities[0]
5146 );
5147 }
5148
5149 #[test]
5150 fn test_kl_counts_returns_temperature_uncertainty() {
5151 let data = u238_single_resonance();
5152 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5153 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5154
5155 let config = UnifiedFitConfig::new(
5156 energies,
5157 vec![data],
5158 vec!["U-238".into()],
5159 350.0,
5160 None,
5161 vec![0.0005],
5162 )
5163 .unwrap()
5164 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5165 .with_fit_temperature(true);
5166
5167 let input = InputData::Counts {
5168 sample_counts: sample,
5169 open_beam_counts: open_beam,
5170 };
5171 let result = fit_spectrum_typed(&input, &config).unwrap();
5172 assert!(result.converged);
5173 let unc = result
5174 .uncertainties
5175 .as_ref()
5176 .expect("KL+temp fit should return density uncertainties");
5177 assert!(
5178 unc[0].is_finite() && unc[0] > 0.0,
5179 "density unc = {}",
5180 unc[0]
5181 );
5182 let t_unc = result
5183 .temperature_k_unc
5184 .expect("KL+temp fit should return temperature uncertainty");
5185 assert!(
5186 t_unc.is_finite() && t_unc > 0.0,
5187 "temperature unc = {t_unc}"
5188 );
5189 }
5190
5191 #[test]
5194 fn test_scale_by_chi2_inflates_joint_poisson_uncertainty() {
5195 let data = u238_single_resonance();
5196 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5197 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5198 let input = InputData::Counts {
5199 sample_counts: sample,
5200 open_beam_counts: open_beam,
5201 };
5202
5203 let make_config = |scale: bool| {
5204 UnifiedFitConfig::new(
5205 energies.clone(),
5206 vec![data.clone()],
5207 vec!["U-238".into()],
5208 350.0,
5209 None,
5210 vec![0.0005],
5211 )
5212 .unwrap()
5213 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5214 .with_fit_temperature(true)
5215 .with_scale_by_chi2(scale)
5216 };
5217
5218 let unscaled = fit_spectrum_typed(&input, &make_config(false)).unwrap();
5219 let scaled = fit_spectrum_typed(&input, &make_config(true)).unwrap();
5220 assert!(unscaled.converged && scaled.converged);
5221
5222 let dpd = unscaled
5225 .deviance_per_dof
5226 .expect("joint-Poisson must report deviance_per_dof");
5227 assert!(dpd.is_finite() && dpd > 0.0, "deviance_per_dof = {dpd}");
5228 let factor = dpd.sqrt();
5229 assert!(
5234 (factor - 1.0).abs() > 0.01,
5235 "expected D/dof to move σ by >1%, got factor {factor}"
5236 );
5237
5238 let t_unscaled = unscaled.temperature_k_unc.expect("σ_T unscaled");
5240 let t_scaled = scaled.temperature_k_unc.expect("σ_T scaled");
5241 let rel_t = (t_scaled - t_unscaled * factor).abs() / (t_unscaled * factor);
5242 assert!(
5243 rel_t < 1e-6,
5244 "σ_T: scaled {t_scaled} must equal unscaled {t_unscaled} × {factor}"
5245 );
5246
5247 let d_unscaled = unscaled.uncertainties.as_ref().expect("density σ unscaled")[0];
5249 let d_scaled = scaled.uncertainties.as_ref().expect("density σ scaled")[0];
5250 let rel_d = (d_scaled - d_unscaled * factor).abs() / (d_unscaled * factor);
5251 assert!(
5252 rel_d < 1e-6,
5253 "density σ: scaled {d_scaled} must equal unscaled {d_unscaled} × {factor}"
5254 );
5255
5256 let default_cfg = UnifiedFitConfig::new(
5259 energies.clone(),
5260 vec![data.clone()],
5261 vec!["U-238".into()],
5262 350.0,
5263 None,
5264 vec![0.0005],
5265 )
5266 .unwrap()
5267 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5268 .with_fit_temperature(true);
5269 let default_run = fit_spectrum_typed(&input, &default_cfg).unwrap();
5270 assert_eq!(
5271 default_run.temperature_k_unc, unscaled.temperature_k_unc,
5272 "default (flag absent) must equal scale_by_chi2=false"
5273 );
5274 }
5275
5276 #[test]
5284 fn test_scale_by_chi2_joint_poisson_underfit_grows_sigma() {
5285 let data = u238_single_resonance();
5286 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5287 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5288 let sample: Vec<f64> = sample
5293 .iter()
5294 .enumerate()
5295 .map(|(i, &c)| {
5296 let k = if i % 2 == 0 { 1.30 } else { 0.70 };
5297 (c * k).round().max(0.0)
5298 })
5299 .collect();
5300 let input = InputData::Counts {
5301 sample_counts: sample,
5302 open_beam_counts: open_beam,
5303 };
5304 let make_config = |scale: bool| {
5305 UnifiedFitConfig::new(
5306 energies.clone(),
5307 vec![data.clone()],
5308 vec!["U-238".into()],
5309 350.0,
5310 None,
5311 vec![0.0005],
5312 )
5313 .unwrap()
5314 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5315 .with_fit_temperature(true)
5316 .with_scale_by_chi2(scale)
5317 };
5318 let unscaled = fit_spectrum_typed(&input, &make_config(false)).unwrap();
5319 let scaled = fit_spectrum_typed(&input, &make_config(true)).unwrap();
5320 assert!(unscaled.converged && scaled.converged);
5321 let dpd = unscaled
5322 .deviance_per_dof
5323 .expect("joint-Poisson must report deviance_per_dof");
5324 assert!(
5325 dpd > 1.0,
5326 "zig-zag counts must under-fit (D/dof > 1), got {dpd}"
5327 );
5328 let t_unscaled = unscaled.temperature_k_unc.expect("σ_T unscaled");
5329 let t_scaled = scaled.temperature_k_unc.expect("σ_T scaled");
5330 assert!(
5331 t_scaled > t_unscaled,
5332 "under-fit (D/dof {dpd} > 1) must GROW σ_T: scaled {t_scaled} \
5333 vs unscaled {t_unscaled}"
5334 );
5335 }
5336
5337 #[test]
5353 fn test_scale_by_chi2_transmission_kl_scales_by_reported_reduced_chi2() {
5354 let data = u238_single_resonance();
5355 let true_density = 0.002;
5356 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5357
5358 let make_config = |scale: bool| {
5359 UnifiedFitConfig::new(
5360 energies.clone(),
5361 vec![data.clone()],
5362 vec!["U-238".into()],
5363 0.0,
5364 None,
5365 vec![0.001],
5366 )
5367 .unwrap()
5368 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5369 .with_scale_by_chi2(scale)
5370 };
5371
5372 let check = |t: Vec<f64>, sigma: Vec<f64>| -> (f64, f64) {
5376 let input = InputData::Transmission {
5377 transmission: t,
5378 uncertainty: sigma,
5379 };
5380 let unscaled = fit_spectrum_typed(&input, &make_config(false)).unwrap();
5381 let scaled = fit_spectrum_typed(&input, &make_config(true)).unwrap();
5382 assert!(unscaled.converged && scaled.converged);
5383 let rcs = scaled.reduced_chi_squared;
5384 assert!(
5385 (rcs - unscaled.reduced_chi_squared).abs() < 1e-9,
5386 "scale flag must not change the reported reduced_chi_squared"
5387 );
5388 assert!(rcs.is_finite() && rcs > 0.0, "reduced_chi_squared = {rcs}");
5389 let factor = rcs.sqrt();
5390 let s_un = unscaled.uncertainties.as_ref().expect("σ unscaled")[0];
5391 let s_sc = scaled.uncertainties.as_ref().expect("σ scaled")[0];
5392 let rel = (s_sc - s_un * factor).abs() / (s_un * factor);
5393 assert!(
5394 rel < 1e-6,
5395 "σ_scaled {s_sc} must equal σ_unscaled {s_un} × √(reduced_chi² {rcs})"
5396 );
5397 (rcs, s_sc / s_un)
5398 };
5399
5400 let zigzag = |k: f64| -> (Vec<f64>, Vec<f64>) {
5408 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
5409 let perturbed: Vec<f64> = t
5410 .iter()
5411 .enumerate()
5412 .map(|(i, &ti)| {
5413 if i % 2 == 0 {
5414 ti * (1.0 + k)
5415 } else {
5416 ti * (1.0 - k)
5417 }
5418 })
5419 .collect();
5420 (perturbed, sigma)
5421 };
5422
5423 let (t_good, s_good) = zigzag(0.003);
5425 let (rcs_good, ratio_good) = check(t_good, s_good);
5426 assert!(
5427 rcs_good < 1.0,
5428 "clean-ish fit should give reduced-χ² < 1, got {rcs_good}"
5429 );
5430 assert!(
5431 ratio_good < 1.0,
5432 "good fit must shrink σ, ratio = {ratio_good}"
5433 );
5434
5435 let (t_poor, s_poor) = zigzag(0.03);
5437 let (rcs_poor, ratio_poor) = check(t_poor, s_poor);
5438 assert!(
5439 rcs_poor > 1.0,
5440 "zig-zag fit should give reduced-χ² > 1, got {rcs_poor}"
5441 );
5442 assert!(
5443 ratio_poor > 1.0,
5444 "poor fit must GROW σ, ratio = {ratio_poor}"
5445 );
5446 }
5447
5448 #[test]
5449 fn test_kl_counts_with_background_returns_uncertainty() {
5450 let data = u238_single_resonance();
5451 let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
5452 let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);
5453
5454 let config = UnifiedFitConfig::new(
5455 energies,
5456 vec![data],
5457 vec!["U-238".into()],
5458 300.0,
5459 None,
5460 vec![0.0005],
5461 )
5462 .unwrap()
5463 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5464 .with_fit_temperature(true)
5465 .with_transmission_background(BackgroundConfig::default());
5466
5467 let input = InputData::Counts {
5468 sample_counts: sample,
5469 open_beam_counts: open_beam,
5470 };
5471 let result = fit_spectrum_typed(&input, &config).unwrap();
5472 assert!(result.converged);
5473 let unc = result
5474 .uncertainties
5475 .as_ref()
5476 .expect("KL+bg fit should return density uncertainties");
5477 assert!(
5478 unc[0].is_finite() && unc[0] > 0.0,
5479 "density unc = {}",
5480 unc[0]
5481 );
5482 }
5483
5484 #[test]
5485 fn test_lm_uncertainty_not_regressed() {
5486 let data = u238_single_resonance();
5487 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
5488 let (t_clean, sigma) = synthetic_transmission(&data, 0.001, &energies);
5489 let t: Vec<f64> = t_clean
5496 .iter()
5497 .enumerate()
5498 .map(|(i, &v)| v * (1.0 + 0.002 * (7.3 * i as f64).sin()))
5499 .collect();
5500
5501 let config = UnifiedFitConfig::new(
5502 energies,
5503 vec![data],
5504 vec!["U-238".into()],
5505 0.0,
5506 None,
5507 vec![0.0005],
5508 )
5509 .unwrap()
5510 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
5511
5512 let input = InputData::Transmission {
5513 transmission: t,
5514 uncertainty: sigma,
5515 };
5516 let result = fit_spectrum_typed(&input, &config).unwrap();
5517 assert!(result.converged);
5518 let unc = result
5519 .uncertainties
5520 .as_ref()
5521 .expect("LM should still return uncertainties");
5522 assert!(unc[0].is_finite() && unc[0] > 0.0);
5523 }
5524
5525 #[test]
5543 fn test_energy_scale_with_temperature_recovers_all_three() {
5544 let data = u238_three_resonances();
5547 let flight_path = 25.0_f64;
5548 let true_density = 0.002;
5549 let true_temp = 340.0;
5550 let true_t0 = 0.6_f64; let true_l_scale = 1.004_f64;
5552 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5554
5555 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5557 &nominal,
5558 true_t0,
5559 true_l_scale,
5560 flight_path,
5561 )
5562 .unwrap();
5563 let (t_obs, sigma) =
5565 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5566
5567 let config = UnifiedFitConfig::new(
5568 nominal,
5569 vec![data],
5570 vec!["U-238".into()],
5571 true_temp - 40.0, None,
5573 vec![true_density],
5574 )
5575 .unwrap()
5576 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5577 .with_fit_temperature(true)
5578 .with_energy_scale(0.0, 1.0, flight_path); let input = InputData::Transmission {
5581 transmission: t_obs,
5582 uncertainty: sigma,
5583 };
5584 let result = fit_spectrum_typed(&input, &config).expect("joint fit runs");
5585 assert!(result.converged, "joint (t0,L_scale,T) fit should converge");
5586
5587 let t0 = result.t0_us.expect("t0_us populated");
5588 let ls = result.l_scale.expect("l_scale populated");
5589 let temp = result.temperature_k.expect("temperature_k populated");
5590 assert!(
5591 (t0 - true_t0).abs() < 0.05,
5592 "t0: fitted={t0}, true={true_t0}"
5593 );
5594 assert!(
5595 (ls - true_l_scale).abs() / true_l_scale < 1e-3,
5596 "l_scale: fitted={ls}, true={true_l_scale}"
5597 );
5598 assert!(
5599 (temp - true_temp).abs() < 3.0,
5600 "temperature: fitted={temp}, true={true_temp}"
5601 );
5602 }
5603
5604 #[test]
5612 fn test_energy_scale_with_temperature_recovers_all_three_kl() {
5613 let data = u238_three_resonances();
5614 let flight_path = 25.0_f64;
5615 let true_density = 0.002;
5616 let true_temp = 340.0;
5617 let true_t0 = 0.6_f64;
5618 let true_l_scale = 1.004_f64;
5619 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5620 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5621 &nominal,
5622 true_t0,
5623 true_l_scale,
5624 flight_path,
5625 )
5626 .unwrap();
5627 let (t_obs, sigma) =
5628 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5629
5630 let config = UnifiedFitConfig::new(
5631 nominal,
5632 vec![data],
5633 vec!["U-238".into()],
5634 true_temp - 40.0,
5635 None,
5636 vec![true_density],
5637 )
5638 .unwrap()
5639 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5640 .with_fit_temperature(true)
5641 .with_energy_scale(0.0, 1.0, flight_path);
5642
5643 let result = fit_spectrum_typed(
5644 &InputData::Transmission {
5645 transmission: t_obs,
5646 uncertainty: sigma,
5647 },
5648 &config,
5649 )
5650 .expect("KL joint fit runs");
5651 assert!(result.converged, "KL joint fit should converge");
5652 let t0 = result.t0_us.expect("t0_us populated");
5653 let ls = result.l_scale.expect("l_scale populated");
5654 let temp = result.temperature_k.expect("temperature_k populated");
5655 assert!(
5656 (t0 - true_t0).abs() < 0.1,
5657 "KL t0: fitted={t0}, true={true_t0}"
5658 );
5659 assert!(
5660 (ls - true_l_scale).abs() / true_l_scale < 2e-3,
5661 "KL l_scale: fitted={ls}, true={true_l_scale}"
5662 );
5663 assert!(
5664 (temp - true_temp).abs() < 5.0,
5665 "KL temperature: fitted={temp}, true={true_temp}"
5666 );
5667 let t_unc = result
5668 .temperature_k_unc
5669 .expect("KL joint fit must report a temperature σ");
5670 assert!(
5671 t_unc.is_finite() && t_unc > 0.0,
5672 "KL joint temperature σ must be finite positive, got {t_unc}"
5673 );
5674 }
5675
5676 #[test]
5682 fn test_energy_scale_with_temperature_recovers_all_three_counts() {
5683 let data = u238_three_resonances();
5684 let flight_path = 25.0_f64;
5685 let true_density = 0.002;
5686 let true_temp = 340.0;
5687 let true_t0 = 0.6_f64;
5688 let true_l_scale = 1.004_f64;
5689 let flux = 1.0e4_f64;
5690 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5691 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5692 &nominal,
5693 true_t0,
5694 true_l_scale,
5695 flight_path,
5696 )
5697 .unwrap();
5698 let (t_true, _) = synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5699 let open_beam = vec![flux; t_true.len()];
5700 let sample_counts: Vec<f64> = t_true.iter().map(|&t| flux * t).collect();
5701
5702 let config = UnifiedFitConfig::new(
5703 nominal,
5704 vec![data],
5705 vec!["U-238".into()],
5706 true_temp - 40.0,
5707 None,
5708 vec![true_density],
5709 )
5710 .unwrap()
5711 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
5712 .with_fit_temperature(true)
5713 .with_energy_scale(0.0, 1.0, flight_path);
5714
5715 let result = fit_spectrum_typed(
5716 &InputData::Counts {
5717 sample_counts,
5718 open_beam_counts: open_beam,
5719 },
5720 &config,
5721 )
5722 .expect("counts joint-Poisson joint fit runs");
5723 assert!(result.converged, "counts joint fit should converge");
5724 let t0 = result.t0_us.expect("t0_us populated");
5725 let ls = result.l_scale.expect("l_scale populated");
5726 let temp = result.temperature_k.expect("temperature_k populated");
5727 assert!(
5728 (t0 - true_t0).abs() < 0.1,
5729 "counts t0: fitted={t0}, true={true_t0}"
5730 );
5731 assert!(
5732 (ls - true_l_scale).abs() / true_l_scale < 2e-3,
5733 "counts l_scale: fitted={ls}, true={true_l_scale}"
5734 );
5735 assert!(
5736 (temp - true_temp).abs() < 5.0,
5737 "counts temperature: fitted={temp}, true={true_temp}"
5738 );
5739 let t_unc = result
5740 .temperature_k_unc
5741 .expect("counts joint fit must report a temperature σ");
5742 assert!(
5743 t_unc.is_finite() && t_unc > 0.0,
5744 "counts joint temperature σ must be finite positive, got {t_unc}"
5745 );
5746 }
5747
5748 #[test]
5766 fn test_energy_scale_temperature_sigma_covariance_plumbing() {
5767 let data = u238_three_resonances();
5768 let flight_path = 25.0_f64;
5769 let true_density = 0.002;
5770 let true_temp = 340.0;
5771 let true_t0 = 0.6_f64;
5772 let true_l_scale = 1.004_f64;
5773 let nominal: Vec<f64> = (0..801).map(|i| 4.0 + (i as f64) * 0.05).collect();
5774 let e_true = nereids_fitting::resolution_calib::corrected_energy_grid(
5775 &nominal,
5776 true_t0,
5777 true_l_scale,
5778 flight_path,
5779 )
5780 .unwrap();
5781 let (t_clean, sigma) =
5782 synthetic_transmission_at_temp(&data, true_density, true_temp, &e_true);
5783 let noise = seeded_gaussian(t_clean.len(), 0x6340_0000_0000_0634);
5786 let t_obs: Vec<f64> = t_clean
5787 .iter()
5788 .zip(sigma.iter())
5789 .zip(noise.iter())
5790 .map(|((&t, &s), &g)| t + s * g)
5791 .collect();
5792
5793 let joint_cfg = UnifiedFitConfig::new(
5797 nominal.clone(),
5798 vec![data.clone()],
5799 vec!["U-238".into()],
5800 true_temp - 40.0,
5801 None,
5802 vec![true_density],
5803 )
5804 .unwrap()
5805 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5806 .with_fit_temperature(true)
5807 .with_energy_scale(0.0, 1.0, flight_path)
5808 .with_tzero_jacobian_method(Some(EnergyScaleJacobianMethod::FiniteDifference));
5809 let joint = fit_spectrum_typed(
5810 &InputData::Transmission {
5811 transmission: t_obs.clone(),
5812 uncertainty: sigma.clone(),
5813 },
5814 &joint_cfg,
5815 )
5816 .unwrap();
5817 let sigma_t_joint = joint.temperature_k_unc.expect("joint σ_T");
5818
5819 let ref_cfg = UnifiedFitConfig::new(
5822 e_true,
5823 vec![data.clone()],
5824 vec!["U-238".into()],
5825 true_temp - 40.0,
5826 None,
5827 vec![true_density],
5828 )
5829 .unwrap()
5830 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
5831 .with_fit_temperature(true);
5832 let reference = fit_spectrum_typed(
5833 &InputData::Transmission {
5834 transmission: t_obs,
5835 uncertainty: sigma.clone(),
5836 },
5837 &ref_cfg,
5838 )
5839 .unwrap();
5840 let sigma_t_ref = reference.temperature_k_unc.expect("reference σ_T");
5841
5842 assert!(sigma_t_joint.is_finite() && sigma_t_joint > 0.0);
5850 assert!(sigma_t_ref.is_finite() && sigma_t_ref > 0.0);
5851 assert!(
5852 sigma_t_joint > sigma_t_ref,
5853 "joint σ_T ({sigma_t_joint:.5}) must strictly exceed no-energy-scale \
5854 σ_T ({sigma_t_ref:.5}): the extra t0/L_scale params are not \
5855 orthogonal to temperature"
5856 );
5857
5858 let p = [
5861 joint.densities[0],
5862 joint.temperature_k.unwrap(),
5863 joint.t0_us.unwrap(),
5864 joint.l_scale.unwrap(),
5865 ];
5866 let model = EnergyScaleTransmissionModel::new(
5867 Arc::new(vec![data]),
5868 Arc::new(vec![0]),
5869 Arc::new(vec![1.0]),
5870 true_temp - 40.0,
5871 nominal.clone(),
5872 flight_path,
5873 2,
5874 3,
5875 None,
5876 )
5877 .with_temperature_index(Some(1))
5878 .expect("distinct temperature index")
5879 .with_jacobian_method(EnergyScaleJacobianMethod::FiniteDifference);
5880
5881 let steps = [1e-4 * p[0].max(1e-6), 1e-4 * p[1].max(1.0), 1e-4, 1e-7];
5886 let n_e = nominal.len();
5887 let mut jac = vec![[0.0f64; 4]; n_e];
5888 for (j, &h) in steps.iter().enumerate() {
5889 let mut pp = p;
5890 let mut pm = p;
5891 pp[j] += h;
5892 pm[j] -= h;
5893 let yp = model.evaluate(&pp).unwrap();
5894 let ym = model.evaluate(&pm).unwrap();
5895 for i in 0..n_e {
5896 jac[i][j] = (yp[i] - ym[i]) / (2.0 * h);
5897 }
5898 }
5899 let mut a = vec![vec![0.0f64; 4]; 4];
5901 for i in 0..n_e {
5902 let w = 1.0 / (sigma[i] * sigma[i]);
5903 for r in 0..4 {
5904 for c in 0..4 {
5905 a[r][c] += jac[i][r] * w * jac[i][c];
5906 }
5907 }
5908 }
5909 let cov = invert_dense(&a).expect("covariance invertible");
5910 let sigma_t_recon = (cov[1][1] * joint.reduced_chi_squared).sqrt();
5912 let rel = (sigma_t_recon - sigma_t_joint).abs() / sigma_t_joint;
5913 assert!(
5914 rel < 1e-3,
5915 "independent covariance reconstruction σ_T={sigma_t_recon:.5} must match \
5916 the solver-reported σ_T={sigma_t_joint:.5} to <0.1%, got {rel:.3e} \
5917 (a mismatch means the solver read T's σ from the wrong covariance slot)"
5918 );
5919 }
5920
5921 #[test]
5924 fn test_spectrum_result_corrected_energies() {
5925 let nominal: Vec<f64> = (0..50).map(|i| 5.0 + i as f64 * 0.3).collect();
5926 let flight_path = 25.0;
5927 let (t0, ls) = (0.4_f64, 1.003_f64);
5928 let kl = nereids_physics::resolution::TOF_FACTOR * flight_path;
5934 let expected: Vec<f64> = nominal
5935 .iter()
5936 .map(|&e| {
5937 let tof = kl / e.sqrt();
5938 (kl * ls / (tof - t0)).powi(2)
5939 })
5940 .collect();
5941
5942 let mk = |t0: Option<f64>, ls: Option<f64>, fp: Option<f64>| SpectrumFitResult {
5943 densities: vec![0.001],
5944 uncertainties: None,
5945 reduced_chi_squared: 1.0,
5946 converged: true,
5947 iterations: 1,
5948 temperature_k: None,
5949 temperature_k_unc: None,
5950 anorm: 1.0,
5951 background: [0.0; 3],
5952 back_d: None,
5953 back_f: None,
5954 t0_us: t0,
5955 l_scale: ls,
5956 energy_scale_flight_path_m: fp,
5957 deviance_per_dof: None,
5958 baseline: None,
5959 baseline_e_ref_ev: None,
5960 warnings: Vec::new(),
5961 };
5962
5963 let got = mk(Some(t0), Some(ls), Some(flight_path))
5966 .corrected_energies(&nominal)
5967 .expect("Some when energy scale fitted")
5968 .unwrap();
5969 assert_eq!(got, expected);
5970 assert!(mk(None, None, None).corrected_energies(&nominal).is_none());
5972 assert!(
5973 mk(Some(t0), None, None)
5974 .corrected_energies(&nominal)
5975 .is_none()
5976 );
5977
5978 assert!(
5984 mk(Some(t0), Some(ls), Some(0.0))
5985 .corrected_energies(&nominal)
5986 .unwrap()
5987 .is_err()
5988 );
5989 assert!(
5990 mk(Some(t0), Some(ls), Some(-25.0))
5991 .corrected_energies(&nominal)
5992 .unwrap()
5993 .is_err()
5994 );
5995 assert!(
5996 mk(Some(t0), Some(-1.0), Some(flight_path))
5997 .corrected_energies(&nominal)
5998 .unwrap()
5999 .is_err()
6000 );
6001 assert!(
6002 mk(Some(t0), Some(f64::NAN), Some(flight_path))
6003 .corrected_energies(&nominal)
6004 .unwrap()
6005 .is_err()
6006 );
6007 assert!(
6011 mk(Some(0.0), Some(1.0), Some(flight_path))
6012 .corrected_energies(&[1.0, f64::NAN, 3.0])
6013 .unwrap()
6014 .is_err()
6015 );
6016 assert!(
6017 mk(Some(t0), Some(ls), Some(flight_path))
6018 .corrected_energies(&[0.0, 1.0])
6019 .unwrap()
6020 .is_err()
6021 );
6022 assert!(
6025 mk(Some(t0), Some(ls), Some(flight_path))
6026 .corrected_energies(&[5.0, 4.0, 6.0])
6027 .unwrap()
6028 .is_err()
6029 );
6030 }
6031
6032 #[test]
6034 fn test_energy_scale_returns_fitted_params() {
6035 let data = u238_single_resonance();
6036 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
6037 let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);
6038
6039 let config = UnifiedFitConfig::new(
6040 energies,
6041 vec![data],
6042 vec!["U-238".into()],
6043 293.6,
6044 None,
6045 vec![0.001],
6046 )
6047 .unwrap()
6048 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
6049 .with_transmission_background(BackgroundConfig::default())
6050 .with_energy_scale(0.0, 1.0, 25.0);
6051
6052 let input = InputData::Transmission {
6053 transmission: t_obs,
6054 uncertainty: sigma,
6055 };
6056 let result = fit_spectrum_typed(&input, &config).unwrap();
6057 assert!(result.converged, "Fit should converge");
6058 assert!(
6059 result.t0_us.is_some(),
6060 "t0_us should be Some when energy-scale is fitted"
6061 );
6062 assert!(
6063 result.l_scale.is_some(),
6064 "l_scale should be Some when energy-scale is fitted"
6065 );
6066 let t0 = result.t0_us.unwrap();
6067 let ls = result.l_scale.unwrap();
6068 assert!(t0.is_finite(), "t0 should be finite, got {t0}");
6070 assert!(ls.is_finite(), "l_scale should be finite, got {ls}");
6071 assert!(t0.abs() < 10.0, "t0 should be within bounds, got {t0}");
6072 assert!(
6073 ls > 0.98 && ls < 1.02,
6074 "l_scale should be within bounds, got {ls}"
6075 );
6076 }
6077
6078 #[test]
6080 fn test_no_energy_scale_returns_none() {
6081 let data = u238_single_resonance();
6082 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
6083 let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);
6084
6085 let config = UnifiedFitConfig::new(
6086 energies,
6087 vec![data],
6088 vec!["U-238".into()],
6089 293.6,
6090 None,
6091 vec![0.001],
6092 )
6093 .unwrap()
6094 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
6095 .with_transmission_background(BackgroundConfig::default());
6096
6097 let input = InputData::Transmission {
6098 transmission: t_obs,
6099 uncertainty: sigma,
6100 };
6101 let result = fit_spectrum_typed(&input, &config).unwrap();
6102 assert!(
6103 result.t0_us.is_none(),
6104 "t0_us should be None without energy-scale"
6105 );
6106 assert!(
6107 result.l_scale.is_none(),
6108 "l_scale should be None without energy-scale"
6109 );
6110 }
6111
6112 #[test]
6122 fn test_joint_poisson_density_recovery_c_5_98() {
6123 let data = u238_single_resonance();
6124 let true_density = 0.0005_f64;
6125 let c = 5.98_f64;
6126 let lam_ob = 1000.0_f64; let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
6128 let (t, _) = synthetic_transmission(&data, true_density, &energies);
6129
6130 let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
6133 let sample_counts: Vec<f64> = t.iter().map(|&ti| c * lam_ob * ti).collect();
6134
6135 let config = UnifiedFitConfig::new(
6136 energies,
6137 vec![data],
6138 vec!["U-238".into()],
6139 0.0,
6140 None,
6141 vec![0.001],
6142 )
6143 .unwrap()
6144 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6145 .with_counts_background(CountsBackgroundConfig {
6146 c,
6147 ..Default::default()
6148 });
6149
6150 let input = InputData::Counts {
6151 sample_counts,
6152 open_beam_counts,
6153 };
6154 let result = fit_spectrum_typed(&input, &config).unwrap();
6155
6156 let d_per_dof = result
6158 .deviance_per_dof
6159 .expect("joint-Poisson solver must populate deviance_per_dof");
6160 assert!(
6161 d_per_dof.is_finite() && d_per_dof >= 0.0,
6162 "deviance_per_dof = {d_per_dof} is not a valid GOF"
6163 );
6164 assert!(
6168 d_per_dof < 0.5,
6169 "noise-free expected-counts fit should give D/dof ≈ 0, got {d_per_dof}"
6170 );
6171 let rel_bias = (result.densities[0] - true_density) / true_density;
6173 assert!(
6174 rel_bias.abs() < 0.05,
6175 "density bias {rel_bias} > 5%: fitted={} truth={true_density}",
6176 result.densities[0]
6177 );
6178 assert!((result.reduced_chi_squared - d_per_dof).abs() < 1e-12);
6180 }
6181
6182 #[test]
6186 fn test_joint_poisson_rejects_alpha_fit() {
6187 let data = u238_single_resonance();
6188 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6189 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
6190 let open_beam_counts: Vec<f64> = vec![500.0; energies.len()];
6191 let sample_counts: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
6192
6193 let config = UnifiedFitConfig::new(
6194 energies,
6195 vec![data],
6196 vec!["U-238".into()],
6197 0.0,
6198 None,
6199 vec![0.001],
6200 )
6201 .unwrap()
6202 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6203 .with_counts_background(CountsBackgroundConfig {
6204 fit_alpha_1: true,
6205 c: 1.0,
6206 ..Default::default()
6207 });
6208
6209 let input = InputData::Counts {
6210 sample_counts,
6211 open_beam_counts,
6212 };
6213 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6214 let msg = err.to_string();
6215 assert!(
6216 msg.contains("fit_alpha_1") || msg.contains("alpha_1"),
6217 "expected alpha_1 rejection message, got: {msg}"
6218 );
6219 }
6220
6221 #[test]
6225 fn test_transmission_poisson_kl_dispatches_to_transmission_path() {
6226 let data = u238_single_resonance();
6227 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6228 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
6229
6230 let config = UnifiedFitConfig::new(
6231 energies,
6232 vec![data],
6233 vec!["U-238".into()],
6234 0.0,
6235 None,
6236 vec![0.001],
6237 )
6238 .unwrap()
6239 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));
6240
6241 let input = InputData::Transmission {
6242 transmission: t,
6243 uncertainty: u,
6244 };
6245 let r = fit_spectrum_typed(&input, &config).unwrap();
6246 assert!(r.deviance_per_dof.is_none());
6250 assert!(r.reduced_chi_squared.is_finite());
6251 }
6252
6253 #[test]
6263 fn test_joint_poisson_with_transmission_background() {
6264 let data = u238_single_resonance();
6265 let true_density = 0.0005_f64;
6266 let true_anorm = 0.85_f64;
6267 let true_ba = 0.03_f64;
6268 let true_bb = -0.01_f64;
6269 let true_bc = 0.0_f64;
6270 let c = 5.98_f64;
6271 let lam_ob = 2000.0_f64;
6272 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
6273 let (t_inner, _) = synthetic_transmission(&data, true_density, &energies);
6274 let t_out: Vec<f64> = t_inner
6275 .iter()
6276 .zip(energies.iter())
6277 .map(|(&ti, &e)| true_anorm * ti + true_ba + true_bb / e.sqrt() + true_bc * e.sqrt())
6278 .collect();
6279 let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
6280 let sample_counts: Vec<f64> = t_out.iter().map(|&ti| c * lam_ob * ti).collect();
6281
6282 let bg = BackgroundConfig {
6283 anorm_init: 1.0,
6284 back_a_init: 0.0,
6285 back_b_init: 0.0,
6286 back_c_init: 0.0,
6287 back_d_init: 0.01,
6288 back_f_init: 1.0,
6289 fit_anorm: true,
6290 fit_back_a: true,
6291 fit_back_b: true,
6292 fit_back_c: true,
6293 fit_back_d: false,
6294 fit_back_f: false,
6295 };
6296
6297 let config = UnifiedFitConfig::new(
6298 energies,
6299 vec![data],
6300 vec!["U-238".into()],
6301 0.0,
6302 None,
6303 vec![0.001],
6304 )
6305 .unwrap()
6306 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6307 .with_counts_background(CountsBackgroundConfig {
6308 c,
6309 ..Default::default()
6310 })
6311 .with_transmission_background(bg)
6312 .with_counts_enable_polish(Some(true));
6320
6321 let input = InputData::Counts {
6322 sample_counts,
6323 open_beam_counts,
6324 };
6325 let r = fit_spectrum_typed(&input, &config).unwrap();
6326
6327 let dpd = r.deviance_per_dof.expect("joint-Poisson must report D/dof");
6337 assert!(
6338 dpd < 1.0,
6339 "D/dof = {dpd} unexpectedly large on noise-free fit — bg params not reaching objective?"
6340 );
6341 assert!(r.densities[0] > 1e-5, "density railed: {}", r.densities[0]);
6343 assert!(
6345 (r.anorm - 1.0).abs() > 0.05,
6346 "A_n did not move from init 1.0 (fitted={})",
6347 r.anorm
6348 );
6349 let bg_moved = r.background.iter().any(|v| v.abs() > 1e-4);
6351 assert!(
6352 bg_moved,
6353 "no bg parameter moved from init 0: {:?}",
6354 r.background
6355 );
6356 }
6357
6358 #[test]
6360 fn test_joint_poisson_requires_back_a_when_back_b_enabled() {
6361 let data = u238_single_resonance();
6362 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6363 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
6364 let ob: Vec<f64> = vec![500.0; energies.len()];
6365 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
6366
6367 let bg = BackgroundConfig {
6368 fit_anorm: true,
6370 fit_back_a: false,
6371 fit_back_b: true,
6372 fit_back_c: false,
6373 fit_back_d: false,
6374 fit_back_f: false,
6375 ..BackgroundConfig::default()
6376 };
6377
6378 let config = UnifiedFitConfig::new(
6379 energies,
6380 vec![data],
6381 vec!["U-238".into()],
6382 0.0,
6383 None,
6384 vec![0.001],
6385 )
6386 .unwrap()
6387 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6388 .with_counts_background(CountsBackgroundConfig {
6389 c: 1.0,
6390 ..Default::default()
6391 })
6392 .with_transmission_background(bg);
6393
6394 let input = InputData::Counts {
6395 sample_counts: s,
6396 open_beam_counts: ob,
6397 };
6398 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6399 let msg = err.to_string();
6400 assert!(
6401 msg.contains("B_A"),
6402 "expected B_A pairing-rule rejection message, got: {msg}"
6403 );
6404 }
6405
6406 #[test]
6411 fn test_joint_poisson_rejects_nonzero_detector_background() {
6412 let data = u238_single_resonance();
6413 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6414 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
6415 let flux: Vec<f64> = vec![500.0; energies.len()];
6416 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
6417 let background: Vec<f64> = vec![5.0; energies.len()];
6418
6419 let config = UnifiedFitConfig::new(
6420 energies,
6421 vec![data],
6422 vec!["U-238".into()],
6423 0.0,
6424 None,
6425 vec![0.001],
6426 )
6427 .unwrap()
6428 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6429 .with_counts_background(CountsBackgroundConfig {
6430 c: 1.0,
6431 ..Default::default()
6432 });
6433
6434 let input = InputData::CountsWithNuisance {
6435 sample_counts: s,
6436 flux,
6437 background,
6438 };
6439 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6440 let msg = err.to_string();
6441 assert!(
6442 msg.contains("B_det"),
6443 "expected deferred-B_det rejection message, got: {msg}"
6444 );
6445 }
6446
6447 #[test]
6449 fn test_joint_poisson_rejects_back_d_f() {
6450 let data = u238_single_resonance();
6451 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6452 let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
6453 let ob: Vec<f64> = vec![500.0; energies.len()];
6454 let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();
6455
6456 let bg = BackgroundConfig {
6457 fit_back_d: true,
6458 ..BackgroundConfig::default()
6459 };
6460
6461 let config = UnifiedFitConfig::new(
6462 energies,
6463 vec![data],
6464 vec!["U-238".into()],
6465 0.0,
6466 None,
6467 vec![0.001],
6468 )
6469 .unwrap()
6470 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6471 .with_counts_background(CountsBackgroundConfig {
6472 c: 1.0,
6473 ..Default::default()
6474 })
6475 .with_transmission_background(bg);
6476
6477 let input = InputData::Counts {
6478 sample_counts: s,
6479 open_beam_counts: ob,
6480 };
6481 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6482 assert!(
6483 err.to_string().contains("BackD"),
6484 "expected BackD/BackF rejection, got: {err}"
6485 );
6486 }
6487
6488 #[test]
6496 fn test_lm_transmission_rejects_partial_back_d_only() {
6497 let data = u238_single_resonance();
6498 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6499 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
6500
6501 let bg = BackgroundConfig {
6502 fit_back_d: true,
6505 fit_back_f: false,
6506 ..BackgroundConfig::default()
6507 };
6508
6509 let config = UnifiedFitConfig::new(
6510 energies,
6511 vec![data],
6512 vec!["U-238".into()],
6513 0.0,
6514 None,
6515 vec![0.001],
6516 )
6517 .unwrap()
6518 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6519 .with_transmission_background(bg);
6520
6521 let input = InputData::Transmission {
6522 transmission: t,
6523 uncertainty: u,
6524 };
6525 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6526 let msg = err.to_string();
6527 assert!(
6528 msg.contains("fit_back_d") && msg.contains("fit_back_f"),
6529 "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
6530 );
6531 }
6532
6533 #[test]
6535 fn test_lm_transmission_rejects_partial_back_f_only() {
6536 let data = u238_single_resonance();
6537 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6538 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
6539
6540 let bg = BackgroundConfig {
6541 fit_back_d: false,
6542 fit_back_f: true,
6543 ..BackgroundConfig::default()
6544 };
6545
6546 let config = UnifiedFitConfig::new(
6547 energies,
6548 vec![data],
6549 vec!["U-238".into()],
6550 0.0,
6551 None,
6552 vec![0.001],
6553 )
6554 .unwrap()
6555 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6556 .with_transmission_background(bg);
6557
6558 let input = InputData::Transmission {
6559 transmission: t,
6560 uncertainty: u,
6561 };
6562 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6563 let msg = err.to_string();
6564 assert!(
6565 msg.contains("fit_back_d") && msg.contains("fit_back_f"),
6566 "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
6567 );
6568 }
6569
6570 #[test]
6572 fn test_transmission_poisson_kl_rejects_partial_back_d_f() {
6573 let data = u238_single_resonance();
6574 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
6575 let (t, u) = synthetic_transmission(&data, 0.0005, &energies);
6576
6577 let bg = BackgroundConfig {
6578 fit_back_d: true,
6579 fit_back_f: false,
6580 ..BackgroundConfig::default()
6581 };
6582
6583 let config = UnifiedFitConfig::new(
6584 energies,
6585 vec![data],
6586 vec!["U-238".into()],
6587 0.0,
6588 None,
6589 vec![0.001],
6590 )
6591 .unwrap()
6592 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6593 .with_transmission_background(bg);
6594
6595 let input = InputData::Transmission {
6596 transmission: t,
6597 uncertainty: u,
6598 };
6599 let err = fit_spectrum_typed(&input, &config).unwrap_err();
6600 assert!(
6601 err.to_string().contains("fit_back_d"),
6602 "expected partial-BackD/F rejection on transmission-PoissonKL path"
6603 );
6604 }
6605
6606 #[test]
6612 fn test_unified_fit_config_fit_energy_range_round_trips() {
6613 let data = u238_single_resonance();
6614 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
6615 let cfg = UnifiedFitConfig::new(
6616 energies,
6617 vec![data],
6618 vec!["U-238".into()],
6619 0.0,
6620 None,
6621 vec![0.001],
6622 )
6623 .unwrap();
6624 assert_eq!(cfg.fit_energy_range(), None);
6625
6626 let cfg = cfg.with_fit_energy_range(Some((5.0, 50.0))).unwrap();
6627 assert_eq!(cfg.fit_energy_range(), Some((5.0, 50.0)));
6628
6629 let cfg = cfg.with_fit_energy_range(None).unwrap();
6631 assert_eq!(cfg.fit_energy_range(), None);
6632 }
6633
6634 #[test]
6638 fn test_unified_fit_config_fit_energy_range_rejects_invalid() {
6639 let data = u238_single_resonance();
6640 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
6641 let cfg = UnifiedFitConfig::new(
6642 energies,
6643 vec![data],
6644 vec!["U-238".into()],
6645 0.0,
6646 None,
6647 vec![0.001],
6648 )
6649 .unwrap();
6650
6651 let err = cfg
6653 .clone()
6654 .with_fit_energy_range(Some((10.0, 5.0)))
6655 .unwrap_err();
6656 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
6657
6658 let err = cfg
6660 .clone()
6661 .with_fit_energy_range(Some((5.0, 5.0)))
6662 .unwrap_err();
6663 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
6664
6665 let err = cfg
6667 .clone()
6668 .with_fit_energy_range(Some((f64::NAN, 5.0)))
6669 .unwrap_err();
6670 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
6671
6672 let err = cfg
6673 .clone()
6674 .with_fit_energy_range(Some((5.0, f64::INFINITY)))
6675 .unwrap_err();
6676 assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
6677 }
6678
6679 #[test]
6684 fn test_fit_energy_range_lm_matches_subset_when_outside_negligible() {
6685 let data = u238_single_resonance();
6686 let true_density = 0.002;
6687 let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
6692 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
6693
6694 let cfg_full = UnifiedFitConfig::new(
6695 energies.clone(),
6696 vec![data.clone()],
6697 vec!["U-238".into()],
6698 0.0,
6699 None,
6700 vec![0.001],
6701 )
6702 .unwrap()
6703 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
6704
6705 let cfg_masked = UnifiedFitConfig::new(
6706 energies.clone(),
6707 vec![data],
6708 vec!["U-238".into()],
6709 0.0,
6710 None,
6711 vec![0.001],
6712 )
6713 .unwrap()
6714 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6715 .with_fit_energy_range(Some((4.0, 8.0)))
6716 .unwrap();
6717
6718 let input = InputData::Transmission {
6719 transmission: t,
6720 uncertainty: sigma,
6721 };
6722 let r_full = fit_spectrum_typed(&input, &cfg_full).unwrap();
6723 let r_masked = fit_spectrum_typed(&input, &cfg_masked).unwrap();
6724 assert!(r_full.converged && r_masked.converged);
6725
6726 let d_full = r_full.densities[0];
6727 let d_masked = r_masked.densities[0];
6728 let rel_err = (d_full - d_masked).abs() / d_full.abs();
6729 assert!(
6730 rel_err < 0.01,
6731 "fit_energy_range LM density {d_masked} should match full-grid {d_full} \
6732 within 1% (got rel_err = {rel_err})"
6733 );
6734 }
6735
6736 #[test]
6743 fn test_fit_energy_range_rejected_on_transmission_poisson_path() {
6744 let data = u238_single_resonance();
6745 let true_density = 0.002;
6746 let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
6747 let (t, sigma) = synthetic_transmission(&data, true_density, &energies);
6748
6749 let cfg = UnifiedFitConfig::new(
6750 energies,
6751 vec![data],
6752 vec!["U-238".into()],
6753 0.0,
6754 None,
6755 vec![0.001],
6756 )
6757 .unwrap()
6758 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6759 .with_fit_energy_range(Some((4.0, 8.0)))
6760 .unwrap();
6761
6762 let input = InputData::Transmission {
6763 transmission: t,
6764 uncertainty: sigma,
6765 };
6766 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
6767 let msg = err.to_string();
6768 assert!(
6769 msg.contains("fit_energy_range")
6770 && (msg.contains("Poisson") || msg.contains("poisson")),
6771 "expected rejection message mentioning fit_energy_range and the \
6772 Poisson-KL path; got: {msg}"
6773 );
6774 }
6775
6776 #[test]
6781 fn test_fit_energy_range_lm_rejects_too_narrow() {
6782 let data = u238_single_resonance();
6783 let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
6787 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
6788 let cfg = UnifiedFitConfig::new(
6789 energies,
6790 vec![data],
6791 vec!["U-238".into()],
6792 0.0,
6793 None,
6794 vec![0.001],
6795 )
6796 .unwrap()
6797 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6798 .with_fit_energy_range(Some((4.6, 4.7)))
6799 .unwrap();
6800 let input = InputData::Transmission {
6801 transmission: t,
6802 uncertainty: sigma,
6803 };
6804 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
6805 let msg = err.to_string();
6806 assert!(
6807 msg.contains("fit_energy_range") && msg.contains("active bin"),
6808 "expected too-narrow-range rejection; got: {msg}"
6809 );
6810 }
6811
6812 #[test]
6815 fn test_fit_energy_range_jp_rejects_too_narrow() {
6816 let data = u238_single_resonance();
6817 let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
6818 let (sample, open_beam) = synthetic_counts(&data, 0.002, &energies, 1000.0);
6819 let cfg = UnifiedFitConfig::new(
6820 energies,
6821 vec![data],
6822 vec!["U-238".into()],
6823 0.0,
6824 None,
6825 vec![0.001],
6826 )
6827 .unwrap()
6828 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6829 .with_fit_energy_range(Some((4.6, 4.7)))
6830 .unwrap();
6831 let input = InputData::Counts {
6832 sample_counts: sample,
6833 open_beam_counts: open_beam,
6834 };
6835 let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
6836 let msg = err.to_string();
6837 assert!(
6838 msg.contains("fit_energy_range") && msg.contains("active bin"),
6839 "expected too-narrow-range rejection; got: {msg}"
6840 );
6841 }
6842
6843 #[test]
6854 fn evaluate_jacobian_and_fisher_gaussian_aux_grid() {
6855 use nereids_physics::resolution::{ResolutionFunction, ResolutionParams};
6856 let data = u238_single_resonance();
6857 let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
6858 let n_e = energies.len();
6859 let config = UnifiedFitConfig::new(
6860 energies,
6861 vec![data],
6862 vec!["U-238".into()],
6863 300.0,
6864 Some(ResolutionFunction::Gaussian(
6865 ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap(),
6866 )),
6867 vec![0.001],
6868 )
6869 .unwrap();
6870 let flux = vec![5000.0; n_e];
6871 let background = vec![10.0; n_e];
6872 let result = evaluate_jacobian_and_fisher(&config, &flux, &background).unwrap();
6873 assert_eq!(result.model_prediction.len(), n_e);
6874 assert!(result.model_prediction.iter().all(|v| v.is_finite()));
6875 assert_eq!(result.param_names.len(), 1, "one free density parameter");
6876 let f00 = result.fisher.get(0, 0);
6879 assert!(f00.is_finite() && f00 > 0.0, "Fisher[0,0] = {f00}");
6880 assert!(result.jacobian.get(0, 0).is_finite());
6881 assert!(result.jacobian.get(n_e - 1, 0).is_finite());
6882 }
6883
6884 #[test]
6894 fn fit_spectrum_typed_energy_scale_lm_recovers_calibration() {
6895 let data = hf178_mlbw_two_resonances(); let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
6897 let true_density = 0.05_f64;
6898 let (t0_true, ls_true) = (1.5_f64, 1.004_f64);
6899 let model = EnergyScaleTransmissionModel::new(
6902 Arc::new(vec![data.clone()]),
6903 Arc::new(vec![0]),
6904 Arc::new(vec![1.0]),
6905 293.6,
6906 energies.clone(),
6907 25.0,
6908 1,
6909 2,
6910 None,
6911 );
6912 let measured = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
6913 let uncertainty = vec![1e-3; energies.len()];
6914 let config = UnifiedFitConfig::new(
6915 energies,
6916 vec![data],
6917 vec!["Hf-178".into()],
6918 293.6,
6919 None,
6920 vec![0.04], )
6922 .unwrap()
6923 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
6924 .with_energy_scale(0.0, 1.0, 25.0); let input = InputData::Transmission {
6926 transmission: measured,
6927 uncertainty,
6928 };
6929 let result = fit_spectrum_typed(&input, &config).unwrap();
6930 let t0 = result
6931 .t0_us
6932 .expect("t0_us populated when fit_energy_scale=true");
6933 let ls = result
6934 .l_scale
6935 .expect("l_scale populated when fit_energy_scale=true");
6936 assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
6937 assert!(result.converged, "energy-scale LM fit should converge");
6938 assert!(
6939 (result.densities[0] - true_density).abs() / true_density < 0.10,
6940 "density: fitted={}, true={true_density}",
6941 result.densities[0]
6942 );
6943 }
6944
6945 #[test]
6950 fn fit_spectrum_typed_energy_scale_counts_kl_seeds_via_proxy() {
6951 let data = hf178_mlbw_two_resonances();
6952 let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
6953 let true_density = 0.05_f64;
6954 let (t0_true, ls_true) = (1.0_f64, 1.003_f64);
6955 let model = EnergyScaleTransmissionModel::new(
6956 Arc::new(vec![data.clone()]),
6957 Arc::new(vec![0]),
6958 Arc::new(vec![1.0]),
6959 293.6,
6960 energies.clone(),
6961 25.0,
6962 1,
6963 2,
6964 None,
6965 );
6966 let t = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
6967 let flux: Vec<f64> = vec![5000.0; energies.len()];
6971 let background: Vec<f64> = vec![0.0; energies.len()];
6972 let sample: Vec<f64> = t
6973 .iter()
6974 .zip(flux.iter())
6975 .map(|(&ti, &fi)| fi * ti)
6976 .collect();
6977 let config = UnifiedFitConfig::new(
6978 energies,
6979 vec![data],
6980 vec!["Hf-178".into()],
6981 293.6,
6982 None,
6983 vec![0.04],
6984 )
6985 .unwrap()
6986 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
6987 .with_energy_scale(0.0, 1.0, 25.0);
6988 let input = InputData::CountsWithNuisance {
6989 sample_counts: sample,
6990 flux,
6991 background,
6992 };
6993 let result = fit_spectrum_typed(&input, &config).unwrap();
6994 let t0 = result
6995 .t0_us
6996 .expect("t0_us populated when fit_energy_scale=true");
6997 let ls = result
6998 .l_scale
6999 .expect("l_scale populated when fit_energy_scale=true");
7000 assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
7001 }
7002
7003 #[test]
7009 fn validate_precomputed_work_cross_sections_error_branches() {
7010 use nereids_physics::transmission::WorkingGridLayout;
7011 let data = u238_single_resonance();
7012 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
7013 let n_e = energies.len();
7014 let n_work = n_e + 2; let good_layout = || {
7016 Arc::new(WorkingGridLayout {
7017 energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
7018 data_indices: (0..n_e).collect(),
7019 })
7020 };
7021 let cfg = |work_xs: Vec<Vec<f64>>, layout: Arc<WorkingGridLayout>| {
7022 UnifiedFitConfig::new(
7023 energies.clone(),
7024 vec![data.clone()],
7025 vec!["U-238".into()],
7026 0.0,
7027 None,
7028 vec![0.001],
7029 )
7030 .unwrap()
7031 .with_precomputed_cross_sections(Arc::new(vec![vec![1.0f64; n_e]])) .with_precomputed_work_cross_sections(Arc::new(work_xs), layout)
7033 };
7034 let expect =
7035 |c: &UnifiedFitConfig, needle: &str| match validate_precomputed_cross_sections(c) {
7036 Err(PipelineError::ShapeMismatch(m)) => {
7037 assert!(m.contains(needle), "expected {needle:?}, got: {m}")
7038 }
7039 other => panic!("expected ShapeMismatch({needle:?}), got {other:?}"),
7040 };
7041 expect(&cfg(vec![], good_layout()), "must not be empty");
7043 expect(
7045 &cfg(vec![vec![1.0; n_work - 1]], good_layout()),
7046 "working-grid energies",
7047 );
7048 let mut nan_row = vec![1.0; n_work];
7050 nan_row[3] = f64::NAN;
7051 expect(&cfg(vec![nan_row], good_layout()), "non-finite");
7052 expect(
7054 &cfg(vec![vec![1.0; n_work], vec![1.0; n_work]], good_layout()),
7055 "rows but",
7056 );
7057 let short = Arc::new(WorkingGridLayout {
7059 energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
7060 data_indices: (0..n_e - 1).collect(),
7061 });
7062 expect(&cfg(vec![vec![1.0; n_work]], short), "layout maps");
7063 let oor = Arc::new(WorkingGridLayout {
7065 energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
7066 data_indices: {
7067 let mut v: Vec<usize> = (0..n_e).collect();
7068 v[0] = n_work + 5;
7069 v
7070 },
7071 });
7072 expect(&cfg(vec![vec![1.0; n_work]], oor), "out of");
7073 }
7074
7075 #[test]
7080 fn evaluate_jacobian_and_fisher_identity_and_temperature_paths() {
7081 let data = u238_single_resonance();
7082 let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
7083 let n_e = energies.len();
7084 let flux = vec![5000.0; n_e];
7085 let background = vec![10.0; n_e];
7086 let cfg_id = UnifiedFitConfig::new(
7088 energies.clone(),
7089 vec![data.clone()],
7090 vec!["U-238".into()],
7091 300.0,
7092 None,
7093 vec![0.001],
7094 )
7095 .unwrap();
7096 let r_id = evaluate_jacobian_and_fisher(&cfg_id, &flux, &background).unwrap();
7097 assert_eq!(r_id.model_prediction.len(), n_e);
7098 let f00 = r_id.fisher.get(0, 0);
7099 assert!(
7100 f00.is_finite() && f00 > 0.0,
7101 "identity-path Fisher[0,0]={f00}"
7102 );
7103 let cfg_t = UnifiedFitConfig::new(
7106 energies,
7107 vec![data],
7108 vec!["U-238".into()],
7109 300.0,
7110 None,
7111 vec![0.001],
7112 )
7113 .unwrap()
7114 .with_fit_temperature(true);
7115 let r_t = evaluate_jacobian_and_fisher(&cfg_t, &flux, &background).unwrap();
7116 assert_eq!(
7117 r_t.param_names.len(),
7118 2,
7119 "density + temperature free params"
7120 );
7121 assert!(r_t.model_prediction.iter().all(|v| v.is_finite()));
7122 }
7123
7124 #[test]
7129 fn fit_spectrum_typed_grouped_precomputed_collapses_by_groups() {
7130 use nereids_core::types::{Isotope, IsotopeGroup};
7131 let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
7132 let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);
7133 let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
7134 let true_density = 0.0005_f64;
7135 let sample = phys_transmission::SampleParams::new(
7136 0.0,
7137 vec![
7138 (rd1.clone(), true_density * 0.6),
7139 (rd2.clone(), true_density * 0.4),
7140 ],
7141 )
7142 .unwrap();
7143 let transmission = phys_transmission::forward_model(&energies, &sample, None).unwrap();
7144 let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();
7145 let per_member = phys_transmission::broadened_cross_sections(
7147 &energies,
7148 &[rd1.clone(), rd2.clone()],
7149 0.0,
7150 None,
7151 None,
7152 )
7153 .unwrap();
7154 let iso1 = Isotope::new(92, 235).unwrap();
7155 let iso2 = Isotope::new(92, 238).unwrap();
7156 let group = IsotopeGroup::custom("U".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();
7157 let config = UnifiedFitConfig::new(
7158 energies,
7159 vec![rd1.clone()],
7160 vec!["placeholder".into()],
7161 0.0,
7162 None,
7163 vec![0.001],
7164 )
7165 .unwrap()
7166 .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
7167 .unwrap()
7168 .with_precomputed_cross_sections(Arc::new(per_member))
7169 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
7170 let input = InputData::Transmission {
7171 transmission,
7172 uncertainty,
7173 };
7174 let result = fit_spectrum_typed(&input, &config).unwrap();
7175 assert_eq!(result.densities.len(), 1, "one group density");
7176 assert!(
7177 (result.densities[0] - true_density).abs() / true_density < 0.01,
7178 "grouped+precomputed density: fitted={}, true={true_density}",
7179 result.densities[0]
7180 );
7181 }
7182
7183 #[test]
7187 fn peak_match_energy_scale_seed_rejects_out_of_bounds() {
7188 let data = hf178_mlbw_two_resonances();
7189 let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
7190 let model = EnergyScaleTransmissionModel::new(
7191 Arc::new(vec![data.clone()]),
7192 Arc::new(vec![0]),
7193 Arc::new(vec![1.0]),
7194 293.6,
7195 energies.clone(),
7196 25.0,
7197 1,
7198 2,
7199 None,
7200 );
7201 let t_obs = model.evaluate(&[0.05, 0.0, 1.03]).unwrap();
7203 let config = UnifiedFitConfig::new(
7204 energies.clone(),
7205 vec![data],
7206 vec!["Hf-178".into()],
7207 293.6,
7208 None,
7209 vec![0.05],
7210 )
7211 .unwrap()
7212 .with_energy_scale(0.0, 1.0, 25.0);
7213 assert!(
7214 peak_match_energy_scale_seed(
7215 &t_obs,
7216 config.energies(),
7217 &config,
7218 25.0,
7219 (-10.0, 10.0),
7220 (0.99, 1.01),
7221 )
7222 .is_none(),
7223 "an out-of-bounds calibration fit must be rejected (cold-start fallback)"
7224 );
7225 }
7226
7227 #[test]
7233 fn test_count_free_params_drops_fixed_densities() {
7234 let data = u238_single_resonance();
7235 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
7236 let base = UnifiedFitConfig::new(
7237 energies,
7238 vec![data.clone(), data],
7239 vec!["a".into(), "b".into()],
7240 300.0,
7241 None,
7242 vec![0.001, 0.001],
7243 )
7244 .unwrap()
7245 .with_fit_temperature(true);
7246
7247 assert_eq!(count_free_params(&base), 3);
7249 assert_eq!(count_free_params(&base.clone().with_fix_densities(true)), 1);
7251 let one_fixed = base.with_density_free(vec![false, true]).unwrap();
7253 assert_eq!(count_free_params(&one_fixed), 2);
7254 }
7255
7256 #[test]
7259 fn test_with_density_free_validates_and_normalises() {
7260 let data = u238_single_resonance();
7261 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
7262 let cfg = UnifiedFitConfig::new(
7263 energies,
7264 vec![data],
7265 vec!["a".into()],
7266 300.0,
7267 None,
7268 vec![0.001],
7269 )
7270 .unwrap();
7271 assert!(matches!(
7273 cfg.clone().with_density_free(vec![true, false]),
7274 Err(FitConfigError::DensityCountMismatch { .. })
7275 ));
7276 let all_free = cfg.with_density_free(vec![true]).unwrap();
7278 assert_eq!(all_free.n_free_density_params(), 1);
7279 assert!(!all_free.density_is_fixed(0));
7280 }
7281
7282 #[test]
7288 fn test_with_groups_rejects_prior_density_freeze() {
7289 use nereids_core::types::{Isotope, IsotopeGroup};
7290 let data = u238_single_resonance();
7291 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
7292 let iso = Isotope::new(92, 238).unwrap();
7293 let group = IsotopeGroup::custom("U-238".into(), vec![(iso, 1.0)]).unwrap();
7294
7295 let cfg = UnifiedFitConfig::new(
7296 energies,
7297 vec![data.clone()],
7298 vec!["placeholder".into()],
7299 300.0,
7300 None,
7301 vec![0.001],
7302 )
7303 .unwrap()
7304 .with_fix_densities(true);
7306 assert!(cfg.density_is_fixed(0), "mask set before grouping");
7307
7308 assert!(matches!(
7310 cfg.with_groups(&[(&group, &[data])], vec![0.001]),
7311 Err(FitConfigError::DensityFreezeBeforeGroups)
7312 ));
7313 }
7314
7315 #[test]
7318 fn test_freeze_after_groups_succeeds() {
7319 use nereids_core::types::{Isotope, IsotopeGroup};
7320 let data = u238_single_resonance();
7321 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.1).collect();
7322 let iso = Isotope::new(92, 238).unwrap();
7323 let group = IsotopeGroup::custom("U-238".into(), vec![(iso, 1.0)]).unwrap();
7324
7325 let cfg = UnifiedFitConfig::new(
7326 energies,
7327 vec![data.clone()],
7328 vec!["placeholder".into()],
7329 300.0,
7330 None,
7331 vec![0.001],
7332 )
7333 .unwrap()
7334 .with_groups(&[(&group, &[data])], vec![0.001])
7335 .unwrap()
7336 .with_fix_densities(true);
7337 assert_eq!(cfg.n_density_params(), 1);
7338 assert!(
7339 cfg.density_is_fixed(0),
7340 "group density frozen when freeze follows grouping"
7341 );
7342 assert_eq!(cfg.n_free_density_params(), 0);
7343 }
7344
7345 #[test]
7349 fn test_all_frozen_no_free_param_rejected() {
7350 let data = u238_single_resonance();
7351 let true_density = 0.0005;
7352 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7353 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, 300.0, &energies);
7354 let config = UnifiedFitConfig::new(
7355 energies,
7356 vec![data],
7357 vec!["U-238".into()],
7358 300.0,
7359 None,
7360 vec![true_density],
7361 )
7362 .unwrap()
7363 .with_fix_densities(true);
7365 assert_eq!(count_free_params(&config), 0);
7366 let input = InputData::Transmission {
7367 transmission: t,
7368 uncertainty: sigma,
7369 };
7370 assert!(matches!(
7371 fit_spectrum_typed(&input, &config),
7372 Err(PipelineError::InvalidParameter(_))
7373 ));
7374 }
7375
7376 #[test]
7381 fn test_fix_densities_recovers_temperature_holds_density() {
7382 let data = u238_single_resonance();
7383 let true_density = 0.0005;
7384 let true_temp = 350.0;
7385 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7386 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
7387
7388 let config = UnifiedFitConfig::new(
7389 energies,
7390 vec![data],
7391 vec!["U-238".into()],
7392 300.0, None,
7394 vec![true_density], )
7396 .unwrap()
7397 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
7398 .with_fit_temperature(true)
7399 .with_fix_densities(true);
7400
7401 let input = InputData::Transmission {
7402 transmission: t,
7403 uncertainty: sigma,
7404 };
7405 let result = fit_spectrum_typed(&input, &config).unwrap();
7406 assert!(result.converged, "T-only fit should converge");
7407
7408 assert_eq!(
7410 result.densities[0], true_density,
7411 "frozen density must not move"
7412 );
7413 let fitted_temp = result
7415 .temperature_k
7416 .expect("temperature_k should be Some when fit_temperature=true");
7417 assert!(
7418 (fitted_temp - true_temp).abs() < 1.0,
7419 "temperature: fitted={fitted_temp}, true={true_temp}"
7420 );
7421
7422 let t_unc = result
7428 .temperature_k_unc
7429 .expect("temperature_k_unc should be Some for a converged T fit");
7430 assert!(
7431 t_unc.is_finite() && t_unc > 0.0,
7432 "frozen-density thermometry must report a finite positive temperature σ, got {t_unc}"
7433 );
7434 let dens_unc = result
7435 .uncertainties
7436 .as_ref()
7437 .expect("converged fit has density uncertainties");
7438 assert!(
7439 dens_unc[0].is_nan(),
7440 "frozen density must report NaN σ (no covariance column), got {}",
7441 dens_unc[0]
7442 );
7443 }
7444
7445 #[test]
7453 fn test_fix_densities_kl_reports_temperature_uncertainty() {
7454 let data = u238_single_resonance();
7455 let true_density = 0.0005;
7456 let true_temp = 350.0;
7457 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7458 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
7459
7460 let config = UnifiedFitConfig::new(
7461 energies,
7462 vec![data],
7463 vec!["U-238".into()],
7464 300.0,
7465 None,
7466 vec![true_density],
7467 )
7468 .unwrap()
7469 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7470 .with_fit_temperature(true)
7471 .with_fix_densities(true);
7472
7473 let input = InputData::Transmission {
7474 transmission: t,
7475 uncertainty: sigma,
7476 };
7477 let result = fit_spectrum_typed(&input, &config).unwrap();
7478 assert!(result.converged, "KL T-only fit should converge");
7479 assert_eq!(
7480 result.densities[0], true_density,
7481 "frozen density must not move"
7482 );
7483 let fitted_temp = result.temperature_k.expect("temperature_k should be Some");
7484 assert!(
7485 (fitted_temp - true_temp).abs() < 5.0,
7486 "KL temperature: fitted={fitted_temp}, true={true_temp}"
7487 );
7488 let t_unc = result
7490 .temperature_k_unc
7491 .expect("KL frozen-density fit must report a temperature σ, not None");
7492 assert!(
7493 t_unc.is_finite() && t_unc > 0.0,
7494 "KL frozen-density temperature σ must be finite positive, got {t_unc}"
7495 );
7496 }
7497
7498 #[test]
7505 fn test_extract_result_maps_uncertainties_past_frozen_density() {
7506 let data = u238_single_resonance();
7507 let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
7508 let config = UnifiedFitConfig::new(
7510 energies,
7511 vec![data.clone(), data],
7512 vec!["a".into(), "b".into()],
7513 300.0,
7514 None,
7515 vec![0.001, 0.002],
7516 )
7517 .unwrap()
7518 .with_fit_temperature(true)
7519 .with_density_free(vec![false, true])
7520 .unwrap();
7521
7522 let result = LmResult {
7525 chi_squared: 1.0,
7526 reduced_chi_squared: 1.0,
7527 iterations: 5,
7528 converged: true,
7529 params: vec![0.001, 0.002, 350.0],
7530 covariance: Some(lm::FlatMatrix::zeros(2, 2)),
7531 uncertainties: Some(vec![0.02, 4.0]),
7532 };
7533
7534 let extracted = extract_result(&config, &result, 2, &[1, 2], None, None).unwrap();
7535 let unc = extracted
7536 .uncertainties
7537 .expect("converged fit surfaces density uncertainties");
7538 assert!(
7539 unc[0].is_nan(),
7540 "frozen density d0 must report NaN σ, got {}",
7541 unc[0]
7542 );
7543 assert_eq!(unc[1], 0.02, "free density d1 σ must come from free slot 0");
7544 assert_eq!(
7545 extracted.temperature_k_unc,
7546 Some(4.0),
7547 "temperature σ must come from free slot 1, not the missing full index 2"
7548 );
7549 }
7550
7551 #[test]
7555 fn test_fix_densities_biased_density_gives_finite_documented_bias() {
7556 let data = u238_single_resonance();
7557 let true_density = 0.0005;
7558 let true_temp = 350.0;
7559 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7560 let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
7561
7562 let config = UnifiedFitConfig::new(
7563 energies,
7564 vec![data],
7565 vec!["U-238".into()],
7566 300.0,
7567 None,
7568 vec![0.9 * true_density], )
7570 .unwrap()
7571 .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
7572 .with_fit_temperature(true)
7573 .with_fix_densities(true);
7574
7575 let input = InputData::Transmission {
7576 transmission: t,
7577 uncertainty: sigma,
7578 };
7579 let result = fit_spectrum_typed(&input, &config).unwrap();
7580 let fitted_temp = result.temperature_k.expect("temperature_k should be Some");
7581 assert_eq!(result.densities[0], 0.9 * true_density);
7583 assert!(
7584 fitted_temp.is_finite() && fitted_temp > 1.0,
7585 "biased-density fit must yield a finite temperature, got {fitted_temp}"
7586 );
7587 }
7588
7589 const BL_TRUE: [f64; 3] = [1.02, -0.03, 0.01];
7594
7595 fn baseline_at(e: f64, e_ref: f64, b: &[f64; 3]) -> f64 {
7596 let z = (e / e_ref).ln();
7597 b[0] + b[1] * z + b[2] * z * z
7598 }
7599
7600 fn apply_truth_baseline(t: &[f64], energies: &[f64]) -> Vec<f64> {
7604 let e_ref = nereids_fitting::transmission_model::baseline_reference_energy(energies);
7605 let out: Vec<f64> = t
7606 .iter()
7607 .zip(energies.iter())
7608 .map(|(&ti, &e)| ti * baseline_at(e, e_ref, &BL_TRUE))
7609 .collect();
7610 let max_rel = t
7611 .iter()
7612 .zip(out.iter())
7613 .map(|(&a, &b)| ((b - a) / a.max(1e-12)).abs())
7614 .fold(0.0_f64, f64::max);
7615 assert!(
7616 max_rel > 0.01,
7617 "non-vacuity pre-check: injected baseline moved the data by only \
7618 {max_rel:.2e} (max relative) — the closed loop would not \
7619 distinguish baseline-on from baseline-off"
7620 );
7621 out
7622 }
7623
7624 #[test]
7625 fn baseline_lm_closed_loop_recovers_temperature_and_coefficients() {
7626 let data = u238_single_resonance();
7627 let true_density = 0.002;
7628 let true_temp = 600.0;
7629 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7630 let (t_pure, _) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
7631 let measured = apply_truth_baseline(&t_pure, &energies);
7632 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
7633 let e_ref = nereids_fitting::transmission_model::baseline_reference_energy(&energies);
7634
7635 let config = UnifiedFitConfig::new(
7639 energies,
7640 vec![data],
7641 vec!["U-238".into()],
7642 500.0,
7643 None,
7644 vec![true_density],
7645 )
7646 .unwrap()
7647 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7648 .with_fit_temperature(true)
7649 .with_fix_densities(true)
7650 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
7651
7652 let input = InputData::Transmission {
7653 transmission: measured,
7654 uncertainty: sigma,
7655 };
7656 let result = fit_spectrum_typed(&input, &config).unwrap();
7657 assert!(result.converged, "baseline LM closed loop should converge");
7658
7659 let fitted_temp = result.temperature_k.expect("temperature fitted");
7660 assert!(
7661 (fitted_temp - true_temp).abs() < 10.0,
7662 "temperature: fitted={fitted_temp}, true={true_temp}"
7663 );
7664
7665 let b = result.baseline.expect("baseline fitted");
7666 for (i, (&fitted, &truth)) in b.iter().zip(BL_TRUE.iter()).enumerate() {
7667 assert!(
7668 (fitted - truth).abs() < 1e-2,
7669 "b{i}: fitted={fitted}, true={truth}"
7670 );
7671 }
7672 let bl_cfg = MultiplicativeBaselineConfig::default();
7673 for (i, (&fitted, bounds)) in b
7674 .iter()
7675 .zip([bl_cfg.b0_bounds, bl_cfg.b1_bounds, bl_cfg.b2_bounds].iter())
7676 .enumerate()
7677 {
7678 assert!(
7679 fitted >= bounds.0 && fitted <= bounds.1,
7680 "b{i} = {fitted} escaped its bounds {bounds:?}"
7681 );
7682 }
7683 let reported_e_ref = result.baseline_e_ref_ev.expect("e_ref reported");
7684 assert!(
7685 (reported_e_ref - e_ref).abs() < 1e-12,
7686 "reported E_ref {reported_e_ref} != geometric midpoint {e_ref}"
7687 );
7688 assert!(
7689 result.warnings.is_empty(),
7690 "no degenerate trio here — warnings should be empty, got {:?}",
7691 result.warnings
7692 );
7693 }
7694
7695 #[test]
7696 fn baseline_kl_transmission_closed_loop() {
7697 let data = u238_single_resonance();
7698 let true_density = 0.002;
7699 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7700 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
7701 let measured = apply_truth_baseline(&t_pure, &energies);
7702 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
7703
7704 let config = UnifiedFitConfig::new(
7705 energies,
7706 vec![data],
7707 vec!["U-238".into()],
7708 0.0,
7709 None,
7710 vec![0.001],
7711 )
7712 .unwrap()
7713 .with_solver(SolverConfig::PoissonKL(PoissonConfig {
7714 max_iter: 500,
7715 ..PoissonConfig::default()
7716 }))
7717 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
7718
7719 let input = InputData::Transmission {
7720 transmission: measured,
7721 uncertainty: sigma,
7722 };
7723 let result = fit_spectrum_typed(&input, &config).unwrap();
7724 assert!(result.converged, "baseline KL closed loop should converge");
7725 let fitted = result.densities[0];
7726 assert!(
7727 (fitted - true_density).abs() / true_density < 0.02,
7728 "density: fitted={fitted}, true={true_density}"
7729 );
7730 let b = result.baseline.expect("baseline fitted");
7731 for (i, (&fit_b, &truth)) in b.iter().zip(BL_TRUE.iter()).enumerate() {
7732 assert!(
7733 (fit_b - truth).abs() < 1e-2,
7734 "b{i}: fitted={fit_b}, true={truth}"
7735 );
7736 }
7737 }
7738
7739 #[test]
7740 fn baseline_counts_jp_closed_loop() {
7741 let data = u238_single_resonance();
7742 let true_density = 0.002;
7743 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7744 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
7745 let t_baselined = apply_truth_baseline(&t_pure, &energies);
7746 let i0 = 10_000.0;
7747 let open_beam: Vec<f64> = vec![i0; t_baselined.len()];
7748 let sample: Vec<f64> = t_baselined
7749 .iter()
7750 .map(|&v| (v * i0).round().max(0.0))
7751 .collect();
7752
7753 let config = UnifiedFitConfig::new(
7754 energies,
7755 vec![data],
7756 vec!["U-238".into()],
7757 0.0,
7758 None,
7759 vec![0.001],
7760 )
7761 .unwrap()
7762 .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
7763 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
7764
7765 let input = InputData::Counts {
7766 sample_counts: sample,
7767 open_beam_counts: open_beam,
7768 };
7769 let result = fit_spectrum_typed(&input, &config).unwrap();
7770 assert!(result.converged, "baseline JP closed loop should converge");
7771 assert!(
7772 result.deviance_per_dof.is_some(),
7773 "counts-KL path reports deviance"
7774 );
7775 let fitted = result.densities[0];
7776 assert!(
7777 (fitted - true_density).abs() / true_density < 0.02,
7778 "density: fitted={fitted}, true={true_density}"
7779 );
7780 let b = result.baseline.expect("baseline fitted");
7781 for (i, (&fit_b, &truth)) in b.iter().zip(BL_TRUE.iter()).enumerate() {
7782 assert!(
7783 (fit_b - truth).abs() < 1e-2,
7784 "b{i}: fitted={fit_b}, true={truth}"
7785 );
7786 }
7787 }
7788
7789 #[test]
7793 fn baseline_rejects_free_anorm_on_all_paths() {
7794 let data = u238_single_resonance();
7795 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
7796 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
7797
7798 let base_config = UnifiedFitConfig::new(
7799 energies.clone(),
7800 vec![data],
7801 vec!["U-238".into()],
7802 0.0,
7803 None,
7804 vec![0.001],
7805 )
7806 .unwrap()
7807 .with_transmission_background(BackgroundConfig::default())
7809 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
7810
7811 let transmission_input = InputData::Transmission {
7812 transmission: t.clone(),
7813 uncertainty: sigma,
7814 };
7815 let counts_input = InputData::Counts {
7816 sample_counts: t.iter().map(|&v| (v * 1000.0).round()).collect(),
7817 open_beam_counts: vec![1000.0; energies.len()],
7818 };
7819
7820 for (label, input, solver) in [
7821 (
7822 "LM transmission",
7823 &transmission_input,
7824 SolverConfig::LevenbergMarquardt(LmConfig::default()),
7825 ),
7826 (
7827 "KL transmission",
7828 &transmission_input,
7829 SolverConfig::PoissonKL(PoissonConfig::default()),
7830 ),
7831 (
7832 "joint-Poisson counts",
7833 &counts_input,
7834 SolverConfig::PoissonKL(PoissonConfig::default()),
7835 ),
7836 ] {
7837 let config = base_config.clone().with_solver(solver);
7838 let err = fit_spectrum_typed(input, &config)
7839 .expect_err(&format!("{label}: free Anorm + baseline must be rejected"));
7840 let msg = err.to_string();
7841 assert!(
7842 msg.contains("Anorm") && msg.contains("fit_anorm = false"),
7843 "{label}: rejection must name the degeneracy and the fix, got: {msg}"
7844 );
7845 }
7846 }
7847
7848 #[test]
7852 fn baseline_with_fixed_anorm_abc_accepted() {
7853 let data = u238_single_resonance();
7854 let true_density = 0.002;
7855 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
7856 let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
7857 let measured = apply_truth_baseline(&t_pure, &energies);
7858 let sigma: Vec<f64> = measured.iter().map(|&v| 0.01 * v.max(0.01)).collect();
7859
7860 let config = UnifiedFitConfig::new(
7861 energies,
7862 vec![data],
7863 vec!["U-238".into()],
7864 0.0,
7865 None,
7866 vec![0.001],
7867 )
7868 .unwrap()
7869 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7870 .with_transmission_background(BackgroundConfig {
7871 fit_anorm: false,
7872 ..BackgroundConfig::default()
7873 })
7874 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
7875
7876 let input = InputData::Transmission {
7877 transmission: measured,
7878 uncertainty: sigma,
7879 };
7880 let result = fit_spectrum_typed(&input, &config).unwrap();
7881 assert!(result.converged, "fixed-Anorm + ABC + baseline should fit");
7882 assert!(
7883 (result.anorm - 1.0).abs() < f64::EPSILON,
7884 "Anorm was frozen at 1.0, got {}",
7885 result.anorm
7886 );
7887 assert!(result.baseline.is_some(), "baseline block reported");
7888 let b = result.baseline.unwrap();
7891 assert!(
7892 (b[0] - BL_TRUE[0]).abs() < 0.02,
7893 "b0: fitted={}, true={}",
7894 b[0],
7895 BL_TRUE[0]
7896 );
7897 }
7898
7899 #[test]
7900 fn baseline_validation_rejects_bad_configs() {
7901 let data = u238_single_resonance();
7902 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
7903 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
7904 let input = InputData::Transmission {
7905 transmission: t,
7906 uncertainty: sigma,
7907 };
7908 let mk_config = |bl: MultiplicativeBaselineConfig| {
7909 UnifiedFitConfig::new(
7910 energies.clone(),
7911 vec![data.clone()],
7912 vec!["U-238".into()],
7913 0.0,
7914 None,
7915 vec![0.001],
7916 )
7917 .unwrap()
7918 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7919 .with_multiplicative_baseline(bl)
7920 };
7921
7922 let err = fit_spectrum_typed(
7924 &input,
7925 &mk_config(MultiplicativeBaselineConfig {
7926 b0_init: f64::NAN,
7927 ..MultiplicativeBaselineConfig::default()
7928 }),
7929 )
7930 .expect_err("NaN b0_init must be rejected");
7931 assert!(err.to_string().contains("b0_init must be finite"));
7932
7933 let err = fit_spectrum_typed(
7935 &input,
7936 &mk_config(MultiplicativeBaselineConfig {
7937 b1_bounds: (0.05, -0.05),
7938 ..MultiplicativeBaselineConfig::default()
7939 }),
7940 )
7941 .expect_err("reversed b1_bounds must be rejected");
7942 assert!(err.to_string().contains("b1_bounds"));
7943
7944 let err = fit_spectrum_typed(
7946 &input,
7947 &mk_config(MultiplicativeBaselineConfig {
7948 b0_init: 1.5,
7949 ..MultiplicativeBaselineConfig::default()
7950 }),
7951 )
7952 .expect_err("out-of-bounds b0_init must be rejected");
7953 assert!(err.to_string().contains("outside"));
7954
7955 let wide_energies: Vec<f64> = (0..51).map(|i| 1.0 * 1.1f64.powi(i)).collect();
7960 let (wt, ws) = synthetic_transmission(&data, 0.002, &wide_energies);
7961 let wide_input = InputData::Transmission {
7962 transmission: wt,
7963 uncertainty: ws,
7964 };
7965 let wide_config = UnifiedFitConfig::new(
7966 wide_energies,
7967 vec![data.clone()],
7968 vec!["U-238".into()],
7969 0.0,
7970 None,
7971 vec![0.001],
7972 )
7973 .unwrap()
7974 .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
7975 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
7976 b0_init: 0.3,
7977 b1_init: 0.4,
7978 b0_bounds: (0.1, 1.1),
7979 b1_bounds: (-1.0, 1.0),
7980 ..MultiplicativeBaselineConfig::default()
7981 });
7982 let err = fit_spectrum_typed(&wide_input, &wide_config)
7983 .expect_err("non-positive initial B(E) must be rejected");
7984 assert!(err.to_string().contains("not strictly"), "got: {err}");
7985 }
7986
7987 #[test]
7988 fn degenerate_trio_produces_warning() {
7989 let data = u238_single_resonance();
7990 let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.2).collect();
7991
7992 let trio = UnifiedFitConfig::new(
7993 energies.clone(),
7994 vec![data.clone()],
7995 vec!["U-238".into()],
7996 300.0,
7997 None,
7998 vec![0.001],
7999 )
8000 .unwrap()
8001 .with_fit_temperature(true)
8002 .with_transmission_background(BackgroundConfig::default());
8003 let w = degenerate_normalization_warning(&trio)
8004 .expect("free Anorm + free T + free density must warn");
8005 assert!(w.contains("degenerate"), "warning names the failure: {w}");
8006
8007 let frozen_density = trio.clone().with_fix_densities(true);
8009 assert!(degenerate_normalization_warning(&frozen_density).is_none());
8010 let no_temp = trio.clone().with_fit_temperature(false);
8011 assert!(degenerate_normalization_warning(&no_temp).is_none());
8012 let fixed_anorm = trio.clone().with_transmission_background(BackgroundConfig {
8013 fit_anorm: false,
8014 ..BackgroundConfig::default()
8015 });
8016 assert!(degenerate_normalization_warning(&fixed_anorm).is_none());
8017
8018 let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
8021 let input = InputData::Transmission {
8022 transmission: t,
8023 uncertainty: sigma,
8024 };
8025 let config = trio.with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
8026 max_iter: 5,
8027 ..LmConfig::default()
8028 }));
8029 let result = fit_spectrum_typed(&input, &config).unwrap();
8030 assert!(
8031 result.warnings.iter().any(|w| w.contains("degenerate")),
8032 "fit result must carry the degenerate-trio warning, got {:?}",
8033 result.warnings
8034 );
8035 }
8036
8037 #[test]
8038 fn count_free_params_includes_baseline_flags() {
8039 let data = u238_single_resonance();
8040 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
8041 let base = UnifiedFitConfig::new(
8042 energies,
8043 vec![data],
8044 vec!["U-238".into()],
8045 0.0,
8046 None,
8047 vec![0.001],
8048 )
8049 .unwrap();
8050 let n0 = count_free_params(&base);
8051
8052 let all_free = base
8053 .clone()
8054 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8055 assert_eq!(count_free_params(&all_free), n0 + 3);
8056
8057 let one_frozen = base
8058 .clone()
8059 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
8060 fit_b1: false,
8061 ..MultiplicativeBaselineConfig::default()
8062 });
8063 assert_eq!(count_free_params(&one_frozen), n0 + 2);
8064
8065 let frozen = base
8066 .clone()
8067 .with_multiplicative_baseline(MultiplicativeBaselineConfig {
8068 fit_b0: false,
8069 fit_b1: false,
8070 fit_b2: false,
8071 ..MultiplicativeBaselineConfig::default()
8072 });
8073 assert_eq!(count_free_params(&frozen), n0);
8074 }
8075
8076 #[test]
8077 fn evaluate_jacobian_and_fisher_rejects_baseline() {
8078 let data = u238_single_resonance();
8079 let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
8080 let n = energies.len();
8081 let config = UnifiedFitConfig::new(
8082 energies,
8083 vec![data],
8084 vec!["U-238".into()],
8085 0.0,
8086 None,
8087 vec![0.001],
8088 )
8089 .unwrap()
8090 .with_multiplicative_baseline(MultiplicativeBaselineConfig::default());
8091 let err = evaluate_jacobian_and_fisher(&config, &vec![1000.0; n], &vec![0.0; n])
8094 .err()
8095 .expect("research Fisher helper must reject a baseline config");
8096 assert!(err.to_string().contains("multiplicative"), "got: {err}");
8097 }
8098
8099 #[test]
8104 fn baseline_init_positivity_scoped_to_fit_window() {
8105 let data = u238_single_resonance();
8106 let energies: Vec<f64> = (0..61)
8110 .map(|i| 1e-3 * 10f64.powf(i as f64 / 10.0))
8111 .collect();
8112 let bl = MultiplicativeBaselineConfig {
8113 b0_init: 0.9,
8114 b2_init: -0.05,
8115 ..MultiplicativeBaselineConfig::default()
8116 };
8117 let base = UnifiedFitConfig::new(
8118 energies,
8119 vec![data],
8120 vec!["U-238".into()],
8121 0.0,
8122 None,
8123 vec![0.001],
8124 )
8125 .unwrap()
8126 .with_multiplicative_baseline(bl);
8127
8128 let err = validate_multiplicative_baseline(&base)
8130 .expect_err("negative B_init at unmasked edge bins must reject");
8131 assert!(err.to_string().contains("not strictly positive"), "{err}");
8132
8133 let windowed = base
8135 .with_fit_energy_range(Some((0.5, 2.0)))
8136 .expect("valid range");
8137 validate_multiplicative_baseline(&windowed)
8138 .expect("inits positive inside the fit window must be accepted");
8139 }
8140
8141 #[test]
8154 fn build_transmission_model_no_temp_matches_forward_model_bit_exact() {
8155 use nereids_physics::transmission::SampleParams;
8156 let data = u238_single_resonance();
8157 let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
8158 let density = 0.002;
8159
8160 let config = UnifiedFitConfig::new(
8162 energies.clone(),
8163 vec![data.clone()],
8164 vec!["U-238".into()],
8165 293.6,
8166 None,
8167 vec![density],
8168 )
8169 .unwrap();
8170 let model = build_transmission_model(&config, 1, None).unwrap();
8171 let t_model = model.evaluate(&[density]).unwrap();
8172 let sample = SampleParams::new(293.6, vec![(data.clone(), density)]).unwrap();
8173 let t_fwd = phys_transmission::forward_model(&energies, &sample, None).unwrap();
8174 assert_eq!(t_model.len(), t_fwd.len());
8175 for (i, (a, b)) in t_model.iter().zip(t_fwd.iter()).enumerate() {
8176 assert_eq!(
8177 a.to_bits(),
8178 b.to_bits(),
8179 "no-resolution arm, bin {i}: {a:e} != {b:e}"
8180 );
8181 }
8182
8183 let res = nereids_physics::resolution::ResolutionFunction::Gaussian(
8185 nereids_physics::resolution::ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap(),
8186 );
8187 let config_res = UnifiedFitConfig::new(
8188 energies.clone(),
8189 vec![data.clone()],
8190 vec!["U-238".into()],
8191 293.6,
8192 Some(res.clone()),
8193 vec![density],
8194 )
8195 .unwrap();
8196 let model_res = build_transmission_model(&config_res, 1, None).unwrap();
8197 let t_model_res = model_res.evaluate(&[density]).unwrap();
8198 let inst = InstrumentParams { resolution: res };
8199 let sample_res = SampleParams::new(293.6, vec![(data.clone(), density)]).unwrap();
8200 let t_fwd_res =
8201 phys_transmission::forward_model(&energies, &sample_res, Some(&inst)).unwrap();
8202 assert_eq!(t_model_res.len(), t_fwd_res.len());
8203 for (i, (a, b)) in t_model_res.iter().zip(t_fwd_res.iter()).enumerate() {
8204 assert_eq!(
8205 a.to_bits(),
8206 b.to_bits(),
8207 "Gaussian-resolution arm, bin {i}: {a:e} != {b:e}"
8208 );
8209 }
8210
8211 let tab_text = "header\n---\n\
8219 5.0 0.0\n\
8220 -0.01 0.0\n\
8221 -0.005 0.5\n\
8222 0.0 1.0\n\
8223 0.005 0.5\n\
8224 0.01 0.0\n\
8225 \n\
8226 200.0 0.0\n\
8227 -0.02 0.0\n\
8228 -0.01 0.5\n\
8229 0.0 1.0\n\
8230 0.01 0.5\n\
8231 0.02 0.0\n";
8232 let tab = nereids_physics::resolution::TabulatedResolution::from_text(tab_text, 25.0)
8233 .expect("synthetic tabulated kernel parses");
8234 let res_tab = nereids_physics::resolution::ResolutionFunction::Tabulated(Arc::new(tab));
8235 let config_tab = UnifiedFitConfig::new(
8236 energies.clone(),
8237 vec![data.clone()],
8238 vec!["U-238".into()],
8239 293.6,
8240 Some(res_tab.clone()),
8241 vec![density],
8242 )
8243 .unwrap();
8244 let model_tab = build_transmission_model(&config_tab, 1, None).unwrap();
8245 let t_model_tab = model_tab.evaluate(&[density]).unwrap();
8246 let inst_tab = InstrumentParams {
8247 resolution: res_tab,
8248 };
8249 let sample_tab = SampleParams::new(293.6, vec![(data, density)]).unwrap();
8250 let t_fwd_tab =
8251 phys_transmission::forward_model(&energies, &sample_tab, Some(&inst_tab)).unwrap();
8252 assert_eq!(t_model_tab.len(), t_fwd_tab.len());
8253 for (i, (a, b)) in t_model_tab.iter().zip(t_fwd_tab.iter()).enumerate() {
8254 assert_eq!(
8255 a.to_bits(),
8256 b.to_bits(),
8257 "tabulated-resolution arm, bin {i}: {a:e} != {b:e}"
8258 );
8259 }
8260 }
8261}