From 1d9874ae07f9ff96080596f2883f9d410ab15b47 Mon Sep 17 00:00:00 2001 From: "Adam J. Jackson" Date: Fri, 31 Jan 2025 13:49:43 +0000 Subject: [PATCH 1/2] Implement alternate _bin_centres_to_bin_edges mode (#363) An optional parameter is provided to change how bin edges are obtained from bin centres: previously the bin edges were constrained to the initial data range, but this can lead to incorrect scaling when performing broadening. Variable-width broadening schemes are now allowed to extrapolate the bin edges in order to get the correct width scaling. Outside of broadening, the default behaviour is unchanged in order to maintain backward compatibility. The changes are tested and some other tweaks made to improve test coverage of these classes. --- CHANGELOG.rst | 14 ++ euphonic/spectra/base.py | 151 +++++++++++++----- .../test/euphonic_test/test_spectrum1d.py | 42 ++++- .../test/euphonic_test/test_spectrum2d.py | 21 +++ 4 files changed, 185 insertions(+), 43 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 16ba8f017..f7ca4f436 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -7,6 +7,20 @@ breaking change: the public classes, functions and type annotations remain importable from euphonic.spectra. +- Bug fixes + + - An optional parameter is provided to change how bin edges are + obtained from bin centres: previously the bin edges were constrained + to the initial data range, but this can lead to incorrect scaling + when performing broadening. Variable-width broadening schemes are + now allowed to extrapolate the bin edges in order to get the correct + width scaling. + + Outside of broadening, the default behaviour is unchanged in order + to maintain backward compatibility. This is likely to be changed + in the next "major version" (i.e. API-breaking release) of + Euphonic. + `v1.4.0.post1 `_ ----------------------------------------------------------------------------------------- diff --git a/euphonic/spectra/base.py b/euphonic/spectra/base.py index fcc8d885a..c08d0b647 100644 --- a/euphonic/spectra/base.py +++ b/euphonic/spectra/base.py @@ -327,11 +327,24 @@ def _bin_edges_to_centres(bin_edges: Quantity) -> Quantity: return bin_edges[:-1] + 0.5*np.diff(bin_edges) @staticmethod - def _bin_centres_to_edges(bin_centres: Quantity) -> Quantity: - return np.concatenate(( - [bin_centres[0]], - (bin_centres[1:] + bin_centres[:-1])/2, - [bin_centres[-1]])) + def _bin_centres_to_edges( + bin_centres: Quantity, + restrict_range: bool = True + ) -> Quantity: + if restrict_range: + return np.concatenate(( + [bin_centres[0]], + (bin_centres[1:] + bin_centres[:-1])/2, + [bin_centres[-1]], + )) + + half_diff = np.diff(bin_centres) * 0.5 + return np.concatenate( + ( + [bin_centres[0] - half_diff[0]], + bin_centres[:-1] + half_diff, + [bin_centres[-1] + half_diff[-1]], + )) @staticmethod def _is_bin_edge(data_length: int, bin_length: int) -> bool: @@ -344,22 +357,29 @@ def _is_bin_edge(data_length: int, bin_length: int) -> bool: f'Unexpected data axis length {data_length} ' f'for bin axis length {bin_length}') - def get_bin_edges(self) -> Quantity: + def get_bin_edges(self, *, restrict_range: bool = True) -> Quantity: """ Get x-axis bin edges. If the size of x_data is one element larger than y_data, x_data is assumed to contain bin edges, but if x_data is the same size, x_data is assumed to contain bin centres and - a conversion is made. In the conversion, the bin edges will - not go outside the existing data bounds so the first and last - bins may be half-size. In addition, each bin edge is assumed - to be halfway between each bin centre, which may not be an - accurate assumption in the case of differently sized bins. + a conversion is made. + + In this case, bin edges are assumed to be halfway between bin centres. + + Parameters + ---------- + restrict_range + If True (default), the bin edges will not go outside the existing + data bounds so the first and last bins may be half-size. This may + be desirable for plotting. Otherwise, the outer bin edges will + extend from the initial data range. """ # Need to use -1 index for y_data so it also works for # Spectrum1DCollection which has y_data shape (n_spectra, bins) if self._is_bin_edge(self.y_data.shape[-1], self.x_data.shape[0]): return self.x_data - return self._bin_centres_to_edges(self.x_data) + return self._bin_centres_to_edges( + self.x_data, restrict_range=restrict_range) def get_bin_centres(self) -> Quantity: """ @@ -376,16 +396,29 @@ def get_bin_centres(self) -> Quantity: return self._bin_edges_to_centres(self.x_data) return self.x_data - def get_bin_widths(self) -> Quantity: + def get_bin_widths(self, *, restrict_range: bool = True) -> Quantity: """ Get x-axis bin widths + + Parameters + ---------- + restrict_range + If True, bin edges will be clamped to the input data range; if + False, they may be extrapolated beyond the initial range of bin + centres. + + False is usually preferable, this default behaviour is for backward + compatibility and will be removed in a future version. + """ - return np.diff(self.get_bin_edges()) + return np.diff(self.get_bin_edges(restrict_range=restrict_range)) def assert_regular_bins(self, message: str = '', rtol: float = 1e-5, - atol: float = 0.) -> None: + atol: float = 0., + restrict_range: bool = True, + ) -> None: """Raise AssertionError if x-axis bins are not evenly spaced. Note that the positional arguments are different from @@ -404,8 +437,16 @@ def assert_regular_bins(self, Absolute tolerance for 'close enough' values. Note this is a bare float and follows the stored units of the bins. + restrict_range + If True, bin edges will be clamped to the input data range; if + False, they may be extrapolated beyond the initial range of bin + centres. + + You should use the value which is consistent with calls to + get_bin_widths() or get_bin_edges(). + """ - bin_widths = self.get_bin_widths() + bin_widths = self.get_bin_widths(restrict_range=restrict_range) # Need to cast to magnitude to use isclose() with atol before Pint 0.21 if not np.all(np.isclose(bin_widths.magnitude, bin_widths.magnitude[0], rtol=rtol, atol=atol)): @@ -678,14 +719,18 @@ def broaden(self: T, x_width, y_broadened = ureg.Quantity(y_broadened, units=self.y_data_unit) elif isinstance(x_width, Callable): - self.assert_regular_bins(message=( - 'Broadening with convolution requires a regular sampling grid.' - )) + self.assert_regular_bins( + message=( + "Broadening with convolution requires a " + "regular sampling grid." + ), + restrict_range=False, + ) y_broadened = variable_width_broadening( - self.get_bin_edges(), + self.get_bin_edges(restrict_range=False), self.get_bin_centres(), x_width, - (self.y_data * self.get_bin_widths()[0]), + (self.y_data * self.get_bin_widths(restrict_range=False)[0]), width_lower_limit=width_lower_limit, width_convention=width_convention, adaptive_error=width_interpolation_error, @@ -933,8 +978,8 @@ def _broaden_spectrum2d_with_function( """ assert axis in ('x', 'y') - bins = spectrum.get_bin_edges(bin_ax=axis) - bin_widths = np.diff(bins) + bins = spectrum.get_bin_edges(bin_ax=axis, restrict_range=False) + bin_widths = spectrum.get_bin_widths(bin_ax=axis, restrict_range=False) if not np.all(np.isclose(bin_widths, bin_widths[0])): bin_width = bin_widths.mean() @@ -979,29 +1024,38 @@ def copy(self: T) -> T: copy.copy(self.x_tick_labels), copy.deepcopy(self.metadata)) - def get_bin_edges(self, bin_ax: Literal['x', 'y'] = 'x') -> Quantity: + def get_bin_edges( + self, + bin_ax: Literal['x', 'y'] = 'x', + *, + restrict_range: bool = True + ) -> Quantity: """ - Get bin edges for the axis specified by bin_ax. If the size of - bin_ax is one element larger than z_data along that axis, bin_ax - is assumed to contain bin edges, but if bin_ax is the same size, - bin_ax is assumed to contain bin centres and a conversion is made. - In the conversion, the bin edges will not go outside the existing - data bounds so the first and last bins may be half-size. In addition, - each bin edge is assumed to be halfway between each bin centre, - which may not be an accurate assumption in the case of differently - sized bins. + Get bin edges for the axis specified by bin_ax. If the size of bin_ax + is one element larger than z_data along that axis, bin_ax is assumed to + contain bin edges. If they are the same size, bin_ax is assumed to + contain bin centres and a conversion is made. + + In this case, bin edges are assumed to be halfway between bin centres. Parameters ---------- bin_ax The axis to get the bin edges for, 'x' or 'y' + + restrict_range + If True (default), the bin edges will not go outside the existing + data bounds so the first and last bins may be half-size. This may + be desirable for plotting. Otherwise, the outer bin edges will + extend from the initial data range. """ enum = {'x': 0, 'y': 1} bin_data = getattr(self, f'{bin_ax}_data') data_ax_len = self.z_data.shape[enum[bin_ax]] if self._is_bin_edge(data_ax_len, bin_data.shape[0]): return bin_data - return self._bin_centres_to_edges(bin_data) + return self._bin_centres_to_edges( + bin_data, restrict_range=restrict_range) def get_bin_centres(self, bin_ax: Literal['x', 'y'] = 'x') -> Quantity: """ @@ -1025,7 +1079,12 @@ def get_bin_centres(self, bin_ax: Literal['x', 'y'] = 'x') -> Quantity: return self._bin_edges_to_centres(bin_data) return bin_data - def get_bin_widths(self, bin_ax: Literal['x', 'y'] = 'x') -> Quantity: + def get_bin_widths( + self, + bin_ax: Literal['x', 'y'] = 'x', + *, + restrict_range: bool = True, + ) -> Quantity: """ Get x-bin widths along specified axis @@ -1033,8 +1092,16 @@ def get_bin_widths(self, bin_ax: Literal['x', 'y'] = 'x') -> Quantity: ---------- bin_ax Axis of interest, 'x' or 'y' + + restrict_range + If True, bin edges will be clamped to the input data range; if + False, they may be extrapolated beyond the initial range of bin + centres. + + False is usually preferable, this default behaviour is for backward + compatibility and will be removed in a future version. """ - bins = self.get_bin_edges(bin_ax) + bins = self.get_bin_edges(bin_ax, restrict_range=restrict_range) return np.diff(bins) def assert_regular_bins( # pylint: disable=arguments-renamed @@ -1042,7 +1109,8 @@ def assert_regular_bins( # pylint: disable=arguments-renamed bin_ax: Literal['x', 'y'], message: str = '', rtol: float = 1e-5, - atol: float = 0. + atol: float = 0., + restrict_range: bool = True, ) -> None: """Raise AssertionError if x-axis bins are not evenly spaced. @@ -1065,8 +1133,15 @@ def assert_regular_bins( # pylint: disable=arguments-renamed Absolute tolerance for 'close enough' values. Note this is a bare float and follows the stored units of the bins. + restrict_range + If True, bin edges will be clamped to the input data range; if + False, they may be extrapolated beyond the initial range of bin + centres. + + You should use the value which is consistent with calls to + get_bin_widths() or get_bin_edges(). """ - bin_widths = self.get_bin_widths(bin_ax) + bin_widths = self.get_bin_widths(bin_ax, restrict_range=restrict_range) if not np.all(np.isclose(bin_widths, bin_widths[0], rtol=rtol, atol=atol)): raise AssertionError( diff --git a/tests_and_analysis/test/euphonic_test/test_spectrum1d.py b/tests_and_analysis/test/euphonic_test/test_spectrum1d.py index 377a83a12..6148b3637 100644 --- a/tests_and_analysis/test/euphonic_test/test_spectrum1d.py +++ b/tests_and_analysis/test/euphonic_test/test_spectrum1d.py @@ -1,5 +1,6 @@ import json +from pint import Quantity import pytest import numpy as np from numpy.polynomial import Polynomial @@ -433,6 +434,11 @@ def test_broaden_uneven_bins_deprecation_raises_value_error(self): with pytest.raises(ValueError): spec1d.broaden(1*ureg('meV')) + with pytest.raises(ValueError, + match="Broadening with convolution requires"): + spec1d.broaden( + lambda energy: np.ones_like(energy) * ureg('meV')) + def test_variable_broadening_consistent(self): """Check variable broadening is consistent with fixed-width method""" spec1d = get_spectrum1d('quartz_666_dos.json') @@ -463,17 +469,43 @@ def fwhm_function(x): check_spectrum1d(variable_broad_sigma, fixed_broad, y_atol=1e-4) check_spectrum1d(variable_broad_sigma, fixed_broad_sigma, y_atol=1e-4) + def test_broaden_consistent_across_bin_definition(self): + """Check broadening results are the same for both bin conventions""" + spec1d_edges = get_spectrum1d('quartz_666_dos.json') + + spec1d_centred = Spectrum1D(spec1d_edges.get_bin_centres(), + spec1d_edges.y_data) + + # Check we really are using both conventions here + assert len(spec1d_edges.x_data) != len(spec1d_centred.x_data) + + fixed_sigma = 1. * ureg("meV") + npt.assert_allclose(spec1d_edges.broaden(x_width=fixed_sigma).y_data, + spec1d_centred.broaden(x_width=fixed_sigma).y_data) + + def sigma_func(energy: Quantity) -> Quantity: + return np.ones_like(energy) * ureg("meV") + npt.assert_allclose(spec1d_edges.broaden(x_width=sigma_func).y_data, + spec1d_centred.broaden(x_width=sigma_func).y_data) + @pytest.mark.parametrize( - 'spectrum1d_file, expected_bin_edges', [ + 'spectrum1d_file, expected_bin_edges, kwargs', [ ('xsq_spectrum1d.json', np.array([0.5, 1., 2., 3., 4., 5., 6.25, 8., 10., 12., 14., 16.5, - 20., 24., 26.])*ureg('1/angstrom')), + 20., 24., 26.]) * ureg('1/angstrom'), + {'restrict_range': True}), + ('xsq_spectrum1d.json', + np.array([0., 1., 2., 3., 4., 5., 6.25, 8., 10., 12., 14., 16.5, + 20., 24., 28.]) * ureg('1/angstrom'), + {'restrict_range': False}), ('xsq_bin_edges_spectrum1d.json', np.array([0., 1., 2., 3., 4., 5., 6., 8., 10., 12., 14., 16., 20., - 24., 28.])*ureg('1/angstrom'))]) - def test_get_bin_edges(self, spectrum1d_file, expected_bin_edges): + 24., 28.]) * ureg('1/angstrom'), + {'restrict_range': True}) + ]) + def test_get_bin_edges(self, spectrum1d_file, expected_bin_edges, kwargs): spec1d = get_spectrum1d(spectrum1d_file) - bin_edges = spec1d.get_bin_edges() + bin_edges = spec1d.get_bin_edges(**kwargs) assert bin_edges.units == expected_bin_edges.units npt.assert_allclose(bin_edges.magnitude, expected_bin_edges.magnitude) diff --git a/tests_and_analysis/test/euphonic_test/test_spectrum2d.py b/tests_and_analysis/test/euphonic_test/test_spectrum2d.py index bf439f753..1dbd90486 100644 --- a/tests_and_analysis/test/euphonic_test/test_spectrum2d.py +++ b/tests_and_analysis/test/euphonic_test/test_spectrum2d.py @@ -1,5 +1,6 @@ import json +from pint import Quantity import pytest import numpy as np from numpy.polynomial import Polynomial @@ -494,6 +495,26 @@ def test_get_bin_centres_with_invalid_data_shape_raises_value_error( with pytest.raises(ValueError): spec2d.get_bin_centres() + def test_broaden_consistent_across_bin_definition(self): + """Check broadening results are the same for both bin conventions""" + spec2d_edges = get_spectrum2d('quartz_bandstructure_sqw.json') + + spec2d_centred = Spectrum2D(spec2d_edges.x_data, + spec2d_edges.get_bin_centres(bin_ax='y'), + spec2d_edges.z_data) + + # Check we really are using both conventions here + assert len(spec2d_edges.y_data) == len(spec2d_centred.y_data) + 1 + + fixed_sigma = 1. * ureg("meV") + npt.assert_allclose(spec2d_edges.broaden(y_width=fixed_sigma).z_data, + spec2d_centred.broaden(y_width=fixed_sigma).z_data) + + def sigma_func(energy: Quantity) -> Quantity: + return np.ones_like(energy) * ureg("meV") + npt.assert_allclose(spec2d_edges.broaden(y_width=sigma_func).z_data, + spec2d_centred.broaden(y_width=sigma_func).z_data) + def test_mul(self): spec = get_spectrum2d('example_spectrum2d.json') From 48f8d114109d5518cd9b3a6c1cef64579ba781eb Mon Sep 17 00:00:00 2001 From: "Adam J. Jackson" Date: Fri, 7 Feb 2025 23:35:59 +0000 Subject: [PATCH 2/2] Create a stub release.yml to enable workflow dispatch (#367) We need _something_ on master to be able to run/test this from other branches. --- .github/workflows/release.yml | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 .github/workflows/release.yml diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml new file mode 100644 index 000000000..09dd4e6ab --- /dev/null +++ b/.github/workflows/release.yml @@ -0,0 +1,15 @@ +name: "Create a release" + +on: + workflow_dispatch: + inputs: + ref: + description: "Target branch" + required: true + type: string + version: + description: "New version number" + required: true + type: string + push: + description: "Push commit, tag and release"