diff --git a/.copier-answers.yml b/.copier-answers.yml index 0653ab1c..867a00e8 100644 --- a/.copier-answers.yml +++ b/.copier-answers.yml @@ -2,7 +2,7 @@ _commit: 8bdcedc _src_path: gh:/EasyScience/EasyProjectTemplate description: A reflectometry python package built on the EasyScience framework. -max_python: '3.12' +max_python: '3.13' min_python: '3.9' orgname: EasyScience packagename: easyreflectometry diff --git a/.github/workflows/python-ci.yml b/.github/workflows/python-ci.yml index cea64fc2..90967cd9 100644 --- a/.github/workflows/python-ci.yml +++ b/.github/workflows/python-ci.yml @@ -28,7 +28,7 @@ jobs: strategy: max-parallel: 4 matrix: - python-version: ['3.11', '3.12'] + python-version: ['3.11', '3.12', '3.13'] os: [ubuntu-latest, macos-latest, windows-2022] runs-on: ${{ matrix.os }} diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index c14850d9..3a4d1036 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,7 +18,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ['3.11','3.12'] + python-version: ['3.11','3.12','3.13'] if: "!contains(github.event.head_commit.message, '[ci skip]')" steps: diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index e4a56d06..078ad7a9 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -102,7 +102,7 @@ Before you submit a pull request, check that it meets these guidelines: 2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring, and add the feature to the list in README.md. -3. The pull request should work for Python, 3.11 and 3.12, and for PyPy. Check +3. The pull request should work for Python, 3.11, 3.12, and 3.13, and for PyPy. Check https://travis-ci.com/easyScience/EasyReflectometryLib/pull_requests and make sure that the tests pass for all supported Python versions. diff --git a/pyproject.toml b/pyproject.toml index 630422ca..2c442dfe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,10 +23,11 @@ classifiers = [ "Programming Language :: Python :: 3 :: Only", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Development Status :: 3 - Alpha" ] -requires-python = ">=3.11,<3.13" +requires-python = ">=3.11,<3.14" dependencies = [ "easyscience @ git+https://github.com/easyscience/corelib.git@develop", @@ -134,11 +135,12 @@ force-single-line = true legacy_tox_ini = """ [tox] isolated_build = True -envlist = py{3.11,3.12} +envlist = py{3.11,3.12,3.13} [gh-actions] python = 3.11: py311 3.12: py312 + 3.13: py313 [gh-actions:env] PLATFORM = ubuntu-latest: linux diff --git a/src/easyreflectometry/fitting.py b/src/easyreflectometry/fitting.py index a33e6570..86df41e3 100644 --- a/src/easyreflectometry/fitting.py +++ b/src/easyreflectometry/fitting.py @@ -31,6 +31,7 @@ def wrapped(*args, **kwargs): self._fit_func = [func_wrapper(m.interface.fit_func, m.unique_name) for m in args] self._models = args self.easy_science_multi_fitter = EasyScienceMultiFitter(args, self._fit_func) + self._fit_results: list[FitResults] | None = None def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup: """ @@ -75,6 +76,7 @@ def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup: dy.append(1 / np.sqrt(variances_masked)) result = self.easy_science_multi_fitter.fit(x, y, weights=dy) + self._fit_results = result new_data = data.copy() for i, _ in enumerate(result): id = refl_nums[i] @@ -99,7 +101,53 @@ def fit_single_data_set_1d(self, data: DataSet1D) -> FitResults: :param data: DataGroup to be fitted to and populated :param method: Optimisation method """ - return self.easy_science_multi_fitter.fit(x=[data.x], y=[data.y], weights=[data.ye])[0] + x_vals = np.asarray(data.x) + y_vals = np.asarray(data.y) + variances = np.asarray(data.ye) + + zero_variance_mask = variances == 0.0 + num_zero_variance = int(np.sum(zero_variance_mask)) + + if num_zero_variance > 0: + warnings.warn( + f'Masked {num_zero_variance} data point(s) in single-dataset fit due to zero variance during fitting.', + UserWarning, + ) + + valid_mask = ~zero_variance_mask + if not np.any(valid_mask): + raise ValueError('Cannot fit single dataset: all points have zero variance.') + + x_vals_masked = x_vals[valid_mask] + y_vals_masked = y_vals[valid_mask] + variances_masked = variances[valid_mask] + + weights = 1.0 / np.sqrt(variances_masked) + result = self.easy_science_multi_fitter.fit(x=[x_vals_masked], y=[y_vals_masked], weights=[weights])[0] + self._fit_results = [result] + return result + + @property + def chi2(self) -> float | None: + """Total chi-squared across all fitted datasets, or None if no fit has been performed.""" + if self._fit_results is None: + return None + return sum(r.chi2 for r in self._fit_results) + + @property + def reduced_chi(self) -> float | None: + """Reduced chi-squared from the most recent fit, or None if no fit has been performed.""" + if self._fit_results is None: + return None + total_chi2 = sum(r.chi2 for r in self._fit_results) + total_points = sum(np.size(r.x) for r in self._fit_results) + n_params = self._fit_results[0].n_pars + total_dof = total_points - n_params + + if total_dof <= 0: + return None + + return total_chi2 / total_dof def switch_minimizer(self, minimizer: AvailableMinimizers) -> None: """ diff --git a/src/easyreflectometry/project.py b/src/easyreflectometry/project.py index 13b4fdda..8c7b6a73 100644 --- a/src/easyreflectometry/project.py +++ b/src/easyreflectometry/project.py @@ -1,5 +1,6 @@ import datetime import json +import logging import os from pathlib import Path from typing import Dict @@ -10,8 +11,8 @@ import numpy as np from easyscience import global_object from easyscience.fitting import AvailableMinimizers -from easyscience.fitting.fitter import DEFAULT_MINIMIZER from easyscience.variable import Parameter +from easyscience.variable.parameter_dependency_resolver import resolve_all_parameter_dependencies from scipp import DataGroup from easyreflectometry.calculators import CalculatorFactory @@ -20,11 +21,9 @@ from easyreflectometry.data.measurement import extract_orso_title from easyreflectometry.data.measurement import load_data_from_orso_file from easyreflectometry.fitting import MultiFitter -from easyreflectometry.model import LinearSpline from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm -from easyreflectometry.model import Pointwise from easyreflectometry.sample import Layer from easyreflectometry.sample import Material from easyreflectometry.sample import MaterialCollection @@ -32,11 +31,13 @@ from easyreflectometry.sample import Sample from easyreflectometry.sample.collections.base_collection import BaseCollection +logger = logging.getLogger(__name__) + Q_MIN = 0.001 Q_MAX = 0.3 Q_RESOLUTION = 500 -DEFAULT_MINIZER = AvailableMinimizers.LMFit_leastsq +DEFAULT_MINIMIZER = AvailableMinimizers.LMFit_leastsq class Project: @@ -48,6 +49,7 @@ def __init__(self): self._calculator = CalculatorFactory() self._experiments: Dict[DataGroup] = {} self._fitter: MultiFitter = None + self._minimizer_selection: AvailableMinimizers = DEFAULT_MINIMIZER self._colors: list[str] = None self._report = None self._q_min: float = None @@ -207,9 +209,8 @@ def models(self, models: ModelCollection) -> None: def fitter(self) -> MultiFitter: if len(self._models): if (self._fitter is None) or (self._fitter_model_index != self._current_model_index): - minimizer = self.minimizer self._fitter = MultiFitter(self._models[self._current_model_index]) - self.minimizer = minimizer + self._fitter.easy_science_multi_fitter.switch_minimizer(self._minimizer_selection) self._fitter_model_index = self._current_model_index return self._fitter @@ -225,10 +226,14 @@ def calculator(self, calculator: str) -> None: def minimizer(self) -> AvailableMinimizers: if self._fitter is not None: return self._fitter.easy_science_multi_fitter.minimizer.enum - return DEFAULT_MINIMIZER + return self._minimizer_selection @minimizer.setter def minimizer(self, minimizer: AvailableMinimizers) -> None: + old_name = getattr(self._minimizer_selection, 'name', str(self._minimizer_selection)) + new_name = getattr(minimizer, 'name', str(minimizer)) + logger.info('Minimizer changed from %s to %s (fitter active: %s)', old_name, new_name, self._fitter is not None) + self._minimizer_selection = minimizer if self._fitter is not None: self._fitter.easy_science_multi_fitter.switch_minimizer(minimizer) @@ -386,21 +391,10 @@ def _apply_resolution_function( ) -> None: """Set the resolution function on *model* based on variance data in *experiment*. - Prefers Pointwise when q-resolution (xe) data is present, otherwise falls - back to LinearSpline when reflectivity error (ye) data is present. - :param experiment: The experiment whose variance data drives the choice. :param model: The model whose resolution function is set. """ - if sum(experiment.xe) != 0: - resolution_function = Pointwise(q_data_points=[experiment.x, experiment.y, experiment.xe]) - model.resolution_function = resolution_function - elif sum(experiment.ye) != 0: - resolution_function = LinearSpline( - q_data_points=experiment.x, - fwhm_values=np.sqrt(experiment.ye), - ) - model.resolution_function = resolution_function + model.resolution_function = PercentageFwhm(5.0) def load_new_experiment(self, path: Union[Path, str]) -> None: new_experiment = load_as_dataset(str(path)) @@ -603,6 +597,8 @@ def as_dict(self, include_materials_not_in_model=False): self._as_dict_add_experiments(project_dict) if self.fitter is not None: project_dict['fitter_minimizer'] = self.fitter.easy_science_multi_fitter.minimizer.name + elif self._minimizer_selection is not None: + project_dict['fitter_minimizer'] = self._minimizer_selection.name if self._calculator is not None: project_dict['calculator'] = self._calculator.current_interface_name if self._colors is not None: @@ -641,7 +637,7 @@ def from_dict(self, project_dict: dict): if 'materials_not_in_model' in keys: self._materials.extend(MaterialCollection.from_dict(project_dict['materials_not_in_model'])) if 'fitter_minimizer' in keys: - self.fitter.easy_science_multi_fitter.switch_minimizer(AvailableMinimizers[project_dict['fitter_minimizer']]) + self.minimizer = AvailableMinimizers[project_dict['fitter_minimizer']] else: self._fitter = None if 'experiments' in keys: @@ -649,6 +645,9 @@ def from_dict(self, project_dict: dict): else: self._experiments = {} + # Resolve any pending parameter dependencies (constraints) after all objects are loaded + resolve_all_parameter_dependencies(self) + def _from_dict_extract_experiments(self, project_dict: dict) -> Dict[int, DataSet1D]: experiments = {} for key in project_dict['experiments'].keys(): diff --git a/tests/summary/test_summary.py b/tests/summary/test_summary.py index 65de9ba3..3bbb186e 100644 --- a/tests/summary/test_summary.py +++ b/tests/summary/test_summary.py @@ -133,7 +133,7 @@ def test_experiments_section(self, project: Project) -> None: assert 'No. of data points' in html assert '408' in html assert 'Resolution function' in html - assert 'Pointwise' in html + assert 'PercentageFwhm' in html def test_experiments_section_percentage_fhwm(self, project: Project) -> None: # When diff --git a/tests/test_fitting.py b/tests/test_fitting.py index fdeba385..446f10b9 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -1,12 +1,15 @@ __author__ = 'github.com/arm61' import os +from unittest.mock import MagicMock +import numpy as np import pytest from easyscience.fitting.minimizers.factory import AvailableMinimizers import easyreflectometry from easyreflectometry.calculators import CalculatorFactory +from easyreflectometry.data import DataSet1D from easyreflectometry.data.measurement import load from easyreflectometry.fitting import MultiFitter from easyreflectometry.model import Model @@ -224,3 +227,152 @@ def test_fitting_with_manual_zero_variance(): assert 'R_0_model' in analysed.keys() assert 'SLD_0' in analysed.keys() assert 'success' in analysed.keys() + + +def test_fit_single_data_set_1d_masks_zero_variance_points(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + captured = {} + mock_result = MagicMock() + mock_result.chi2 = 1.0 + mock_result.n_pars = 1 + + def _fake_fit(*, x, y, weights): + captured['x'] = x + captured['y'] = y + captured['weights'] = weights + return [mock_result] + + fitter.easy_science_multi_fitter = MagicMock() + fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit) + + data = DataSet1D( + name='single_dataset', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.01, 0.0, 0.04]), + ) + + with pytest.warns(UserWarning, match='Masked 1 data point\(s\) in single-dataset fit'): + result = fitter.fit_single_data_set_1d(data) + + assert result is mock_result + assert np.allclose(captured['x'][0], np.array([0.01, 0.03])) + assert np.allclose(captured['y'][0], np.array([1.0, 0.6])) + assert np.allclose(captured['weights'][0], np.array([10.0, 5.0])) + + +def test_reduced_chi_uses_global_dof_across_fit_results(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + fit_result_1 = MagicMock() + fit_result_1.chi2 = 10.0 + fit_result_1.x = np.arange(5) + fit_result_1.n_pars = 3 + + fit_result_2 = MagicMock() + fit_result_2.chi2 = 14.0 + fit_result_2.x = np.arange(7) + fit_result_2.n_pars = 3 + + fitter._fit_results = [fit_result_1, fit_result_2] + + expected = (10.0 + 14.0) / ((5 + 7) - 3) + assert fitter.reduced_chi == pytest.approx(expected) + + +def test_fit_single_data_set_1d_all_zero_variance_raises(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + data = DataSet1D( + name='all_zero', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.0, 0.0, 0.0]), + ) + + with pytest.raises(ValueError, match='all points have zero variance'): + fitter.fit_single_data_set_1d(data) + + +def test_chi2_returns_none_before_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + assert fitter.chi2 is None + + +def test_chi2_returns_total_after_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + r1 = MagicMock() + r1.chi2 = 5.0 + r2 = MagicMock() + r2.chi2 = 3.0 + + fitter._fit_results = [r1, r2] + assert fitter.chi2 == pytest.approx(8.0) + + +def test_reduced_chi_returns_none_before_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + assert fitter.reduced_chi is None + + +def test_reduced_chi_returns_none_when_dof_zero(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + r1 = MagicMock() + r1.chi2 = 5.0 + r1.x = np.arange(3) + r1.n_pars = 3 # total_points == n_params => dof == 0 + + fitter._fit_results = [r1] + assert fitter.reduced_chi is None + + +def test_fit_single_data_set_1d_no_zero_variance(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + captured = {} + mock_result = MagicMock() + mock_result.chi2 = 2.0 + mock_result.n_pars = 1 + + def _fake_fit(*, x, y, weights): + captured['x'] = x + captured['y'] = y + captured['weights'] = weights + return [mock_result] + + fitter.easy_science_multi_fitter = MagicMock() + fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit) + + data = DataSet1D( + name='no_zero', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.01, 0.04, 0.09]), + ) + + result = fitter.fit_single_data_set_1d(data) + + assert result is mock_result + assert np.allclose(captured['x'][0], np.array([0.01, 0.02, 0.03])) + assert np.allclose(captured['y'][0], np.array([1.0, 0.8, 0.6])) diff --git a/tests/test_ort_file.py b/tests/test_ort_file.py index fe40e748..c547b1f5 100644 --- a/tests/test_ort_file.py +++ b/tests/test_ort_file.py @@ -56,6 +56,7 @@ def fit_model(load_data): reflectivity = data['data']['R_0'].values scale_factor = 1 / np.max(reflectivity) data['data']['R_0'].values *= scale_factor + data['data']['R_0'].variances *= scale_factor**2 # Create a model for the sample @@ -81,22 +82,30 @@ def fit_model(load_data): sample=multi_sample, scale=1, background=0.000001, - resolution_function=PercentageFwhm(0), + resolution_function=PercentageFwhm(5), name='Multilayer Model', ) # Set the fitting parameters - sio2_layer.roughness.bounds = (3, 12) - sio2_layer.material.sld.bounds = (3.47, 5) - sio2_layer.thickness.bounds = (10, 30) - - subphase.material.sld.bounds = (6, 6.35) - dlipids_layer.thickness.bounds = (30, 60) - dlipids_layer.roughness.bounds = (3, 10) - dlipids_layer.material.sld.bounds = (4, 6) - multi_layer_model.scale.bounds = (0.8, 1.2) - multi_layer_model.background.bounds = (1e-6, 1e-3) + sio2_layer.roughness.min = 3 + sio2_layer.roughness.max = 12 + sio2_layer.material.sld.min = 3.47 + sio2_layer.material.sld.max = 5 + sio2_layer.thickness.min = 10 + sio2_layer.thickness.max = 30 + + subphase.material.sld.min = 6 + dlipids_layer.thickness.min = 30 + dlipids_layer.thickness.max = 60 + dlipids_layer.roughness.min = 3 + dlipids_layer.roughness.max = 10 + dlipids_layer.material.sld.min = 4 + dlipids_layer.material.sld.max = 6 + multi_layer_model.scale.min = 0.8 + multi_layer_model.scale.max = 1.2 + multi_layer_model.background.min = 1e-6 + multi_layer_model.background.max = 1e-3 sio2_layer.roughness.free = True sio2_layer.material.sld.free = True @@ -176,4 +185,4 @@ def test_analyze_reduced_data__fit_model_success(fit_model): def test_analyze_reduced_data__fit_model_reasonable(fit_model): - assert fit_model['reduced_chi'] < 0.01 + assert fit_model['reduced_chi'] < 6.0 diff --git a/tests/test_project.py b/tests/test_project.py index 40febe90..b93df068 100644 --- a/tests/test_project.py +++ b/tests/test_project.py @@ -16,7 +16,6 @@ from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm -from easyreflectometry.model import Pointwise from easyreflectometry.project import Project from easyreflectometry.sample import Layer from easyreflectometry.sample import Material @@ -350,6 +349,7 @@ def test_as_dict(self): keys.sort() assert keys == [ 'calculator', + 'fitter_minimizer', 'info', 'models', 'with_experiments', @@ -594,7 +594,7 @@ def test_load_experiment(self): assert isinstance(project.experiments[5], DataSet1D) assert project.experiments[5].name == 'Example data file from refnx docs' assert project.experiments[5].model == model_5 - assert isinstance(project.models[5].resolution_function, Pointwise) + assert isinstance(project.models[5].resolution_function, PercentageFwhm) assert isinstance(project.models[4].resolution_function, PercentageFwhm) def test_load_experiment_sets_resolution_function_pointwise_when_xe_present(self, tmp_path): @@ -610,10 +610,10 @@ def test_load_experiment_sets_resolution_function_pointwise_when_xe_present(self # Then project.load_experiment_for_model_at_index(str(fpath)) - # Expect Pointwise because xe (q-resolution) is present - from easyreflectometry.model.resolution_functions import Pointwise + # Resolution is always set to PercentageFwhm + from easyreflectometry.model.resolution_functions import PercentageFwhm - assert isinstance(project.models[0].resolution_function, Pointwise) + assert isinstance(project.models[0].resolution_function, PercentageFwhm) def test_load_experiment_sets_linearspline_when_only_ye_present(self, tmp_path): # When @@ -628,10 +628,10 @@ def test_load_experiment_sets_linearspline_when_only_ye_present(self, tmp_path): # Then project.load_experiment_for_model_at_index(str(fpath)) - # Expect LinearSpline because only ye (fwhm from ye) is available - from easyreflectometry.model.resolution_functions import LinearSpline + # Resolution is always set to PercentageFwhm + from easyreflectometry.model.resolution_functions import PercentageFwhm - assert isinstance(project.models[0].resolution_function, LinearSpline) + assert isinstance(project.models[0].resolution_function, PercentageFwhm) def test_experimental_data_at_index(self): # When