diff --git a/docs/docs/user-guide/analysis-workflow/analysis.md b/docs/docs/user-guide/analysis-workflow/analysis.md index 4de074cbe..1eae23d6d 100644 --- a/docs/docs/user-guide/analysis-workflow/analysis.md +++ b/docs/docs/user-guide/analysis-workflow/analysis.md @@ -265,6 +265,9 @@ Now, you can inspect the fitted parameters to see how they have changed during the refinement process, select more parameters to be refined, and perform additional fits as needed. +The equations and interpretation of the reported fit-quality values are +listed on the [Reliability Factors](reliability-factors.md) page. + To plot the measured and calculated data after the fit, you can use the `pattern` method of the `display` object: diff --git a/docs/docs/user-guide/analysis-workflow/reliability-factors.md b/docs/docs/user-guide/analysis-workflow/reliability-factors.md new file mode 100644 index 000000000..8261ef6df --- /dev/null +++ b/docs/docs/user-guide/analysis-workflow/reliability-factors.md @@ -0,0 +1,156 @@ +--- +title: Reliability Factors +icon: material/chart-bell-curve-cumulative +--- + +# :material-chart-bell-curve-cumulative: Reliability Factors + +EasyDiffraction reports several complementary measures of agreement +between observed and calculated diffraction data. They are ratios +internally; values shown with a percent sign are multiplied by 100. + +For the equations below, $y_i^{\mathrm{obs}}$ and $y_i^{\mathrm{calc}}$ +are the observed and calculated intensities, $\sigma_i$ is the standard +uncertainty of the observed intensity, and + +$$ +w_i = \frac{1}{\sigma_i^2} +$$ + +is its inverse-variance weight. The sum runs over the $N$ data points +included in the reported value. + +## R-factor (Rf) + +The unweighted profile R-factor is the absolute difference between +observed and calculated intensities, normalized by the total absolute +observed intensity: + +$$ +R_f = +\frac{\sum_i \left|y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right|} + {\sum_i \left|y_i^{\mathrm{obs}}\right|}. +$$ + +Lower values indicate closer agreement. This factor does not use the +measurement uncertainties, so every data point contributes according to +the magnitude of its absolute residual. + +## Squared-residual R-factor (Rf²) + +The value labelled `Rf²` in the fit summary is the unweighted +root-squared-residual ratio: + +$$ +R_{f^2} = +\left[ +\frac{\sum_i \left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2} + {\sum_i \left(y_i^{\mathrm{obs}}\right)^2} +\right]^{1/2}. +$$ + +Despite the compact `Rf²` label, this is not the algebraic square of +$R_f$. The label indicates that squared intensities and residuals are +used before taking the square root. + +## Weighted R-factor (wR) + +The weighted R-factor is the root ratio of weighted squared residuals: + +$$ +wR = +\left[ +\frac{\sum_i w_i + \left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2} + {\sum_i w_i \left(y_i^{\mathrm{obs}}\right)^2} +\right]^{1/2}, +\qquad +w_i = \frac{1}{\sigma_i^2}. +$$ + +Consequently, a point with a smaller standard uncertainty has more +influence than a less precise point. EasyDiffraction expects standard +uncertainties as input and converts them to inverse-variance weights; it +does not use $\sigma_i$ itself as the weight. + +## Chi-square and reduced chi-square + +The uncertainty-weighted sum of squared residuals is + +$$ +\chi^2 = +\sum_i \left( +\frac{y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}}{\sigma_i} +\right)^2 += \sum_i w_i +\left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2. +$$ + +If $p$ free parameters were fitted, the number of degrees of freedom is +$\nu=N-p$, and the reported goodness-of-fit is the reduced chi-square: + +$$ +\chi_\nu^2 = \frac{\chi^2}{\nu}. +$$ + +A value near 1 means that the size of the residuals is consistent with +the stated standard uncertainties. A much larger value can indicate a +poor model or underestimated uncertainties; a much smaller value can +indicate overestimated uncertainties or an over-flexible model. + +For a joint fit, EasyDiffraction also multiplies each experiment's +squared normalized residuals by its normalized joint-fit weight. Those +experiment weights are normalized so that their sum equals the number of +experiments. + +## Expected weighted profile R-factor + +For powder fits, the expected weighted profile R-factor is + +$$ +wR_{\mathrm{expected}} = +\left[ +\frac{\nu} + {\sum_i w_i \left(y_i^{\mathrm{obs}}\right)^2} +\right]^{1/2}. +$$ + +It is the weighted profile R-factor expected when $\chi_\nu^2=1$. +Therefore, $wR / wR_{\mathrm{expected}} = \sqrt{\chi_\nu^2}$ when the +same data points and weights are used for both values. + +## Bragg R-factor (BR) + +When observed and calculated structure-factor magnitudes are available, +EasyDiffraction can report the Bragg R-factor: + +$$ +BR = +\frac{\sum_h \left|F_h^{\mathrm{obs}}-F_h^{\mathrm{calc}}\right|} + {\sum_h F_h^{\mathrm{obs}}}. +$$ + +Here $h$ indexes reflections and $F_h$ is a structure-factor magnitude. +Lower values indicate closer agreement between observed and calculated +reflection amplitudes. + +## Names and data subsets + +The fit summary uses the short labels `Rf`, `Rf²`, `wR`, and `BR`. Saved +deterministic fit results also expose IUCr-style names: + +| Saved result | Definition and scope | +| -------------------- | ---------------------------------------------------------------- | +| `R_factor_all` | $R_f$ for all included observations | +| `wR_factor_all` | $wR$ for all included observations | +| `R_factor_gt` | $R_f$ for observations satisfying $y_i^{\mathrm{obs}}>3\sigma_i$ | +| `wR_factor_gt` | $wR$ for observations satisfying $y_i^{\mathrm{obs}}>3\sigma_i$ | +| `prof_R_factor` | $R_f$ for all included powder-profile points | +| `prof_wR_factor` | $wR$ for all included powder-profile points | +| `prof_wR_expected` | $wR_{\mathrm{expected}}$ for all included powder-profile points | +| `reduced_chi_square` | $\chi_\nu^2$ for the fitted residual vector | + +Only finite observations with finite calculated values and positive, +finite standard uncertainties are included in saved deterministic +statistics. A metric is unavailable when its denominator is zero or when +it does not apply to the fitted data. diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml index f53fc69fe..3b7ae5379 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -222,6 +222,7 @@ nav: - Structure: user-guide/analysis-workflow/model.md - Experiment: user-guide/analysis-workflow/experiment.md - Analysis: user-guide/analysis-workflow/analysis.md + - Reliability Factors: user-guide/analysis-workflow/reliability-factors.md - Report: user-guide/analysis-workflow/report.md - Tutorials: - Tutorials: tutorials/index.md @@ -235,7 +236,7 @@ nav: - HS pd-neut-cwl: tutorials/refine-hs-hrpt.ipynb - Si pd-neut-tof: tutorials/refine-si-sepd.ipynb - PbSO4 pd-xray-cwl: tutorials/refine-pbso4-xray.ipynb - - LaM(7)O3 P02.1 pd-xray-cwl: tutorials/refine-lam7o3-p021.ipynb + - LaM(7)O3 pd-xray-cwl: tutorials/refine-lam7o3-p021.ipynb - LMO pd-neut-cwl: tutorials/refine-lmo-echidna.ipynb - Without Measured Data: - LBCO pd-neut-cwl: tutorials/simulate-lbco-cwl.ipynb diff --git a/src/easydiffraction/analysis/analysis.py b/src/easydiffraction/analysis/analysis.py index da6f174a6..2334bdd4d 100644 --- a/src/easydiffraction/analysis/analysis.py +++ b/src/easydiffraction/analysis/analysis.py @@ -44,6 +44,7 @@ from easydiffraction.analysis.fit_helpers.bayesian import PosteriorSamples from easydiffraction.analysis.fit_helpers.bayesian import posterior_predictive_cache_key from easydiffraction.analysis.fit_helpers.metrics import calculate_r_factor +from easydiffraction.analysis.fit_helpers.metrics import calculate_weighted_r_factor from easydiffraction.analysis.fit_helpers.reporting import FitResults from easydiffraction.analysis.fitting import Fitter from easydiffraction.analysis.fitting import FitterFitOptions @@ -2035,12 +2036,7 @@ def _weighted_r_factor_or_none( """Return a weighted R factor when inputs are available.""" if observed.size == 0: return None - weights = 1.0 / uncertainties**2 - denominator = float(np.sum(weights * observed**2)) - if denominator <= 0.0: - return None - numerator = float(np.sum(weights * (observed - calculated) ** 2)) - value = np.sqrt(numerator / denominator) + value = calculate_weighted_r_factor(observed, calculated, uncertainties) return float(value) if np.isfinite(value) else None @staticmethod diff --git a/src/easydiffraction/analysis/fit_helpers/metrics.py b/src/easydiffraction/analysis/fit_helpers/metrics.py index af05ba648..c8604cd8e 100644 --- a/src/easydiffraction/analysis/fit_helpers/metrics.py +++ b/src/easydiffraction/analysis/fit_helpers/metrics.py @@ -44,7 +44,7 @@ def calculate_r_factor( def calculate_weighted_r_factor( y_obs: np.ndarray, y_calc: np.ndarray, - weights: np.ndarray, + standard_uncertainties: np.ndarray, ) -> float: """ Calculate weighted R-factor between observed and calculated data. @@ -55,8 +55,9 @@ def calculate_weighted_r_factor( Observed data points. y_calc : np.ndarray Calculated data points. - weights : np.ndarray - Weights for each data point. + standard_uncertainties : np.ndarray + Standard uncertainties for the observed data points. The + inverse-variance weights are calculated as ``1 / sigma**2``. Returns ------- @@ -65,7 +66,8 @@ def calculate_weighted_r_factor( """ y_obs = np.asarray(y_obs) y_calc = np.asarray(y_calc) - weights = np.asarray(weights) + standard_uncertainties = np.asarray(standard_uncertainties) + weights = 1.0 / standard_uncertainties**2 numerator = np.sum(weights * (y_obs - y_calc) ** 2) denominator = np.sum(weights * y_obs**2) return np.sqrt(numerator / denominator) if denominator != 0 else np.nan diff --git a/tests/integration/fitting/test_icsd_structure_cif_import.py b/tests/integration/fitting/test_icsd_structure_cif_import.py new file mode 100644 index 000000000..986a6f2b6 --- /dev/null +++ b/tests/integration/fitting/test_icsd_structure_cif_import.py @@ -0,0 +1,151 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause +"""Regression coverage for importing an underscore-style ICSD CIF.""" + +import numpy as np +import pytest + +from easydiffraction.analysis.calculators.cryspy import CryspyCalculator +from easydiffraction.datablocks.structure.item.factory import StructureFactory + +ZRW2O8_CIF = """\ +data_83267-ICSD +_database_code_ICSD 83267 +_audit_creation_date 1998-06-26 +_chemical_name_systematic +'Zirconium bis(tungstate)' +_chemical_formula_structural +'Zr (W O4)2' +_chemical_formula_sum +'O8 W2 Zr1' +_publ_section_title +'O8 W2 Zr1' +loop_ +_citation_id +_citation_journal_abbrev +_citation_year +_citation_journal_volume +_citation_page_first +_citation_page_last +_citation_journal_id_ASTM +primary 'Science' 1996 272 90 92 SCIEAS +loop_ +_publ_author_name +Mary, T.A.;Evans, J.S.O.;Vogt, T.;Sleight, A.W. +_cell_length_a 9.15993(5) +_cell_length_b 9.15993(5) +_cell_length_c 9.15993(5) +_cell_angle_alpha 90. +_cell_angle_beta 90. +_cell_angle_gamma 90. +_cell_volume 768.56 +_cell_formula_units_Z 4 +_symmetry_space_group_name_H-M 'P 21 3' +_symmetry_Int_Tables_number 198 +_refine_ls_R_factor_all 0.024000 +loop_ +_symmetry_equiv_pos_site_id +_symmetry_equiv_pos_as_xyz + 1 '-z+1/2, -x, y+1/2' + 2 '-y+1/2, -z, x+1/2' + 3 '-x+1/2, -y, z+1/2' + 4 '-z, x+1/2, -y+1/2' + 5 '-y, z+1/2, -x+1/2' + 6 '-x, y+1/2, -z+1/2' + 7 'z+1/2, -x+1/2, -y' + 8 'y+1/2, -z+1/2, -x' + 9 'x+1/2, -y+1/2, -z' + 10 'z, x, y' + 11 'y, z, x' + 12 'x, y, z' +loop_ +_atom_type_symbol +_atom_type_oxidation_number +Zr4+ 4 +W6+ 6 +O2- -2 +loop_ +_atom_site_label +_atom_site_type_symbol +_atom_site_symmetry_multiplicity +_atom_site_Wyckoff_symbol +_atom_site_fract_x +_atom_site_fract_y +_atom_site_fract_z +_atom_site_B_iso_or_equiv +_atom_site_occupancy +Zr1 Zr4+ 4 a 0.0003(4) 0.0003(4) 0.0003(4) 0.010(1) 1. +W1 W6+ 4 a 0.3412(3) 0.3412(3) 0.3412(3) 0.012(1) 1. +W2 W6+ 4 a 0.6008(3) 0.6008(3) 0.6008(3) 0.010(1) 1. +O1 O2- 12 b 0.2071(3) 0.4378(4) 0.4470(3) 0.022(1) 1. +O2 O2- 12 b 0.7876(3) 0.5694(4) 0.5565(3) 0.020(1) 1. +O3 O2- 4 a 0.4916(5) 0.4916(5) 0.4916(5) 0.023(1) 1. +O4 O2- 4 a 0.2336(3) 0.2336(3) 0.2336(3) 0.037(1) 1. +""" + + +def test_icsd_cif_import_preserves_structure_and_ionic_symbols_for_cryspy(): + """Import the ICSD structure and preserve ionic atom types through Cryspy.""" + from cryspy.H_functions_global.function_1_cryspy_objects import str_to_globaln + + structure = StructureFactory.from_cif_str(ZRW2O8_CIF) + + assert structure.name == '83267-icsd' + assert structure.as_cif.startswith('data_83267-icsd\n') + assert structure.space_group.name_h_m.value == 'P 21 3' + np.testing.assert_allclose( + [ + structure.cell.length_a.value, + structure.cell.length_b.value, + structure.cell.length_c.value, + structure.cell.angle_alpha.value, + structure.cell.angle_beta.value, + structure.cell.angle_gamma.value, + ], + [9.15993, 9.15993, 9.15993, 90.0, 90.0, 90.0], + ) + np.testing.assert_allclose( + [ + structure.cell.length_a.uncertainty, + structure.cell.length_b.uncertainty, + structure.cell.length_c.uncertainty, + ], + [0.00005, 0.00005, 0.00005], + ) + + expected_sites = { + 'Zr1': ('Zr4+', 4, 'a', (0.0003, 0.0003, 0.0003), (0.0004,) * 3, 0.010), + 'W1': ('W6+', 4, 'a', (0.3412, 0.3412, 0.3412), (0.0003,) * 3, 0.012), + 'W2': ('W6+', 4, 'a', (0.6008, 0.6008, 0.6008), (0.0003,) * 3, 0.010), + 'O1': ('O2-', 12, 'b', (0.2071, 0.4378, 0.4470), (0.0003, 0.0004, 0.0003), 0.022), + 'O2': ('O2-', 12, 'b', (0.7876, 0.5694, 0.5565), (0.0003, 0.0004, 0.0003), 0.020), + 'O3': ('O2-', 4, 'a', (0.4916, 0.4916, 0.4916), (0.0005,) * 3, 0.023), + 'O4': ('O2-', 4, 'a', (0.2336, 0.2336, 0.2336), (0.0003,) * 3, 0.037), + } + + assert structure.atom_sites.names == list(expected_sites) + for label, expected in expected_sites.items(): + type_symbol, multiplicity, wyckoff, coordinates, coordinate_sus, adp_iso = expected + site = structure.atom_sites[label] + assert site.type_symbol.value == type_symbol + assert site.multiplicity.value == multiplicity + assert site.wyckoff_letter.value == wyckoff + assert site.occupancy.value == 1.0 + np.testing.assert_allclose( + [site.fract_x.value, site.fract_y.value, site.fract_z.value], + coordinates, + ) + np.testing.assert_allclose( + [site.fract_x.uncertainty, site.fract_y.uncertainty, site.fract_z.uncertainty], + coordinate_sus, + ) + assert site.adp_iso.value == pytest.approx(adp_iso) + assert site.adp_iso.uncertainty == pytest.approx(0.001) + + cryspy_cif = CryspyCalculator()._convert_structure_to_cryspy_cif(structure) + cryspy_structure = str_to_globaln(cryspy_cif).items[0] + + assert cryspy_structure.data_name == '83267-icsd' + assert [site.type_symbol for site in cryspy_structure.atom_site.items] == [ + expected[0] for expected in expected_sites.values() + ] diff --git a/tests/unit/easydiffraction/analysis/fit_helpers/test_metrics.py b/tests/unit/easydiffraction/analysis/fit_helpers/test_metrics.py index d87f24614..f696b2988 100644 --- a/tests/unit/easydiffraction/analysis/fit_helpers/test_metrics.py +++ b/tests/unit/easydiffraction/analysis/fit_helpers/test_metrics.py @@ -11,12 +11,12 @@ def test_calculate_r_metrics_and_chi_square(): y_obs = np.array([1.0, 2.0, 3.0]) y_calc = np.array([1.1, 1.9, 2.8]) - weights = np.array([1.0, 2.0, 3.0]) + standard_uncertainties = np.array([1.0, 2.0, 3.0]) residuals = y_obs - y_calc r = M.calculate_r_factor(y_obs, y_calc) rb = M.calculate_rb_factor(y_obs, y_calc) - rw = M.calculate_weighted_r_factor(y_obs, y_calc, weights) + rw = M.calculate_weighted_r_factor(y_obs, y_calc, standard_uncertainties) r2 = M.calculate_r_factor_squared(y_obs, y_calc) chi2 = M.calculate_reduced_chi_square(residuals, num_parameters=1) @@ -28,6 +28,23 @@ def test_calculate_r_metrics_and_chi_square(): assert np.isfinite(chi2) +def test_calculate_weighted_r_factor_uses_inverse_variance_weights(): + from easydiffraction.analysis.fit_helpers import metrics as M + + y_obs = np.array([10.0, 20.0]) + y_calc = np.array([9.0, 16.0]) + standard_uncertainties = np.array([1.0, 2.0]) + + result = M.calculate_weighted_r_factor( + y_obs, + y_calc, + standard_uncertainties, + ) + + # sqrt((1 * 1^2 + 1/4 * 4^2) / (1 * 10^2 + 1/4 * 20^2)) + assert np.isclose(result, np.sqrt(5.0 / 200.0)) + + def test_get_reliability_inputs_collects_arrays_with_default_su(): from easydiffraction.analysis.fit_helpers import metrics as M diff --git a/tests/unit/easydiffraction/analysis/test_analysis_coverage.py b/tests/unit/easydiffraction/analysis/test_analysis_coverage.py index 1c05f4aad..34042ba40 100644 --- a/tests/unit/easydiffraction/analysis/test_analysis_coverage.py +++ b/tests/unit/easydiffraction/analysis/test_analysis_coverage.py @@ -576,10 +576,12 @@ def test_weighted_r_factor_or_none_empty_and_zero_denominator(self): def test_weighted_r_factor_or_none_computes_value(self): from easydiffraction.analysis.analysis import Analysis - observed = np.asarray([10.0, 10.0], dtype=float) - calculated = np.asarray([10.0, 10.0], dtype=float) - uncertainties = np.asarray([1.0, 1.0], dtype=float) - assert Analysis._weighted_r_factor_or_none(observed, calculated, uncertainties) == 0.0 + observed = np.asarray([10.0, 20.0], dtype=float) + calculated = np.asarray([9.0, 16.0], dtype=float) + uncertainties = np.asarray([1.0, 2.0], dtype=float) + value = Analysis._weighted_r_factor_or_none(observed, calculated, uncertainties) + assert value is not None + assert np.isclose(value, np.sqrt(5.0 / 200.0)) def test_expected_weighted_r_factor_guards(self): from easydiffraction.analysis.analysis import Analysis