diff --git a/.prettierignore b/.prettierignore index ec19b5ef1..2c912ce56 100644 --- a/.prettierignore +++ b/.prettierignore @@ -28,6 +28,9 @@ node_modules # Tox .tox +# Git worktrees checked out inside the repository +.worktrees/ + # Misc .benchmarks .cache diff --git a/docs/docs/tutorials/delta_lorentz.ipynb b/docs/docs/tutorials/delta_lorentz.ipynb index 47cf6bd90..128e7434e 100644 --- a/docs/docs/tutorials/delta_lorentz.ipynb +++ b/docs/docs/tutorials/delta_lorentz.ipynb @@ -65,7 +65,7 @@ "id": "0aee03b1", "metadata": {}, "source": [ - "Both `A_0` and `lorentzian_width` are here allowed to vary with Q. We here change a few of them just to show how this impacts the model. The `# noqa` comment is because we are accessing private members of the model (ones beginning with `_`), which is generally discouraged. Because of these changes, in the figure below, the delta function at Q=1.25 Å^-1 is much larger than the other ones, and the Lorentzian at Q=1.75 Å^-1 is much narrower and taller than the other ones." + "Both `A_0` and `lorentzian_width` are here allowed to vary with Q. We here change a few of them just to show how this impacts the model. The `# ruff: ignore[private-member-access]` comment is because we are accessing private members of the model (ones beginning with `_`), which is generally discouraged. Because of these changes, in the figure below, the delta function at Q=1.25 Å^-1 is much larger than the other ones, and the Lorentzian at Q=1.75 Å^-1 is much narrower and taller than the other ones." ] }, { diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml index 3db04ace0..6575f0c95 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -4,7 +4,7 @@ site_url: https://easyscience.github.io/dynamics-lib # Repository repo_url: https://github.com/easyscience/dynamics-lib -edit_uri: edit/develop/docs/ +edit_uri: edit/develop/docs/docs/ # Copyright copyright: © 2025-2026 EasyDynamics diff --git a/pixi.lock b/pixi.lock index f07f02cf8..772d51ce2 100644 --- a/pixi.lock +++ b/pixi.lock @@ -7536,15 +7536,19 @@ packages: requires_dist: - darkdetect - easyscience>=2.5.1 + - h5py - ipykernel - ipympl - ipython - ipywidgets - jupyterlab - matplotlib + - numpy - pixi-kernel - plopp - pooch + - scipp + - scipy - sympy - build ; extra == 'dev' - copier ; extra == 'dev' diff --git a/pixi.toml b/pixi.toml index db46523e1..007abbce4 100644 --- a/pixi.toml +++ b/pixi.toml @@ -119,7 +119,10 @@ docstring-lint-check = 'pydoclint --quiet src/' notebook-lint-check = 'nbqa ruff docs/docs/tutorials/' py-lint-check = 'ruff check src/ tests/ docs/docs/tutorials/' py-format-check = 'ruff format --check src/ tests/ docs/docs/tutorials/' -nonpy-format-check = 'npx prettier --list-different --config=prettierrc.toml --ignore-unknown .' +# Refreshes prettier first, so the local version can never drift behind the one CI installs. +nonpy-format-check = { cmd = 'npx prettier --list-different --config=prettierrc.toml --ignore-unknown .', depends-on = [ + 'prettier-install', +] } nonpy-format-check-modified = 'python tools/nonpy_prettier_modified.py' check = 'pre-commit run --hook-stage manual --all-files' @@ -134,7 +137,9 @@ notebook-lint-fix = 'nbqa ruff --fix docs/docs/tutorials/' py-lint-fix = 'ruff check --fix src/ tests/ docs/docs/tutorials/' py-lint-fix-unsafe = 'ruff check --fix --unsafe-fixes src/ tests/ docs/docs/tutorials/' py-format-fix = 'ruff format src/ tests/ docs/docs/tutorials/' -nonpy-format-fix = 'npx prettier --write --list-different --config=prettierrc.toml --ignore-unknown .' +nonpy-format-fix = { cmd = 'npx prettier --write --list-different --config=prettierrc.toml --ignore-unknown .', depends-on = [ + 'prettier-install', +] } nonpy-format-fix-modified = 'python tools/nonpy_prettier_modified.py --write' success-message = 'echo "✅ All auto-formatting steps completed successfully!"' @@ -178,9 +183,11 @@ cov = { depends-on = [ ######################## notebook-convert = 'jupytext docs/docs/tutorials/*.py --from py:percent --to ipynb' -notebook-strip = 'nbstripout docs/docs/tutorials/*.ipynb' +notebook-strip = 'nbstripout docs/docs/tutorials/**/*.ipynb' notebook-tweak = 'python tools/tweak_notebooks.py docs/docs/tutorials/' -notebook-exec = 'python -m pytest --nbmake docs/docs/tutorials/ --nbmake-timeout=1200 --overwrite --color=yes -n auto -v' +notebook-exec = { cmd = 'python -m pytest --nbmake docs/docs/tutorials/ --nbmake-timeout=1200 --overwrite --color=yes -n auto -v', depends-on = [ + 'prefetch-tutorial-data', +] } notebook-prepare = { depends-on = [ #'notebook-convert', @@ -270,7 +277,10 @@ default-build = 'python -m build' dist-build = 'python -m build --wheel --outdir dist' npm-config = 'npm config set registry https://registry.npmjs.org/' -prettier-install = 'npm install --no-save --no-audit --no-fund prettier prettier-plugin-toml' +# --prefix . keeps the install inside this repository: without it, npm walks up the directory +# tree and a stray package.json in a parent directory can silently pin an old prettier. +# @latest keeps local runs on the same version CI installs. +prettier-install = 'npm install --prefix . --no-save --no-audit --no-fund prettier@latest prettier-plugin-toml@latest' clean-pycache = "find . -type d -name '__pycache__' -prune -exec rm -rf '{}' +" diff --git a/pyproject.toml b/pyproject.toml index 3330304ea..c8ff0397c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,10 @@ classifiers = [ requires-python = '>=3.12' dependencies = [ 'easyscience>=2.5.1', # The base library of the EasyScience framework. 2.5.1 adds fitting.Sampler + 'numpy', # Numerical arrays (used directly throughout the library) + 'scipy', # Numerical routines (convolution, interpolation, special functions) + 'scipp', # Labelled multi-dimensional arrays; backs Experiment data handling + 'h5py', # HDF5 backend for scipp's HDF5 I/O (Experiment.load_hdf5) 'matplotlib', # Plotting (posterior trace, corner, and predictive plots) 'pooch', # Data downloader 'darkdetect', # Detecting dark mode (system-level) @@ -167,7 +171,10 @@ fail_under = 0 # Minimum coverage percentage to pass [tool.pytest.ini_options] addopts = '--import-mode=importlib' -markers = ['fast: mark test as fast (should be run on every push)'] +markers = [ + 'fast: mark test as fast (should be run on every push)', + 'network: mark test as downloading data files (deselect with -m "not network" when offline)', +] testpaths = ['tests'] ######################## diff --git a/src/easydynamics/__init__.py b/src/easydynamics/__init__.py index 0b8e27cbd..de7391225 100644 --- a/src/easydynamics/__init__.py +++ b/src/easydynamics/__init__.py @@ -9,16 +9,16 @@ """ from easydynamics.analysis import Analysis +from easydynamics.analysis import Analysis1d from easydynamics.analysis import BoundsSuggestion from easydynamics.analysis import BoundsSuggestions +from easydynamics.analysis import FitBinding from easydynamics.analysis import MultiQPosteriorSampler from easydynamics.analysis import ParameterAnalysis from easydynamics.analysis import ParameterLabels from easydynamics.analysis import ParameterPosterior from easydynamics.analysis import PosteriorSampler from easydynamics.analysis import PosteriorSummary -from easydynamics.analysis.analysis1d import Analysis1d -from easydynamics.analysis.fit_binding import FitBinding from easydynamics.base_classes import EasyDynamicsBase from easydynamics.base_classes import EasyDynamicsModelBase from easydynamics.convolution import Convolution @@ -42,11 +42,11 @@ from easydynamics.settings import ConvolutionSettings from easydynamics.settings import DetailedBalanceSettings from easydynamics.utils import detailed_balance_factor +from easydynamics.utils import hbar from easydynamics.utils import plot_corner from easydynamics.utils import plot_posterior_predictive from easydynamics.utils import plot_trace from easydynamics.utils import slicerplot_with_residuals -from easydynamics.utils.utils import hbar __all__ = [ 'Analysis', diff --git a/src/easydynamics/analysis/__init__.py b/src/easydynamics/analysis/__init__.py index c6eb02a92..2cb853cef 100644 --- a/src/easydynamics/analysis/__init__.py +++ b/src/easydynamics/analysis/__init__.py @@ -2,6 +2,8 @@ # SPDX-License-Identifier: BSD-3-Clause from easydynamics.analysis.analysis import Analysis +from easydynamics.analysis.analysis1d import Analysis1d +from easydynamics.analysis.fit_binding import FitBinding from easydynamics.analysis.parameter_analysis import ParameterAnalysis from easydynamics.analysis.posterior import BoundsSuggestion from easydynamics.analysis.posterior import BoundsSuggestions @@ -13,8 +15,10 @@ __all__ = [ 'Analysis', + 'Analysis1d', 'BoundsSuggestion', 'BoundsSuggestions', + 'FitBinding', 'MultiQPosteriorSampler', 'ParameterAnalysis', 'ParameterLabels', diff --git a/src/easydynamics/analysis/analysis.py b/src/easydynamics/analysis/analysis.py index 89af3fc87..f06a818bd 100644 --- a/src/easydynamics/analysis/analysis.py +++ b/src/easydynamics/analysis/analysis.py @@ -1,10 +1,12 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import warnings from copy import copy from typing import Any import numpy as np +import plopp as pp import scipp as sc from easyscience.fitting.minimizers.utils import FitResults from easyscience.fitting.multi_fitter import MultiFitter @@ -298,6 +300,10 @@ def rebin( self.instrument_model.clear_Q(confirm=True) self._analysis_list_is_dirty = True + self._owner_index = None + # The cached MultiFitter holds the old Analysis1d objects, and the Sampler binds its data + # at construction, so both are stale after a rebin. + self._invalidate_fitter() def calculate( self, @@ -448,11 +454,6 @@ def plot_data_and_model( self._verify_bool(add_background, 'add_background') self._verify_bool(plot_residuals, 'plot_residuals') - if energy is None: - energy = self.energy - - import plopp as pp - data_and_model = self.data_and_model_to_datagroup( energy=energy, add_background=add_background, @@ -464,7 +465,8 @@ def plot_data_and_model( plot_kwargs_defaults['keep'] = 'energy' plot_kwargs_defaults.update(kwargs) - if plot_residuals: + # Residuals may have been omitted (with a warning) for a custom energy grid. + if plot_residuals and 'Residuals' in data_and_model: fig = slicerplot_with_residuals( data_and_model, residuals_key='Residuals', @@ -531,8 +533,20 @@ def data_and_model_to_datagroup( self._verify_bool(include_components, 'include_components') self._verify_bool(include_residuals, 'include_residuals') + custom_energy = energy is not None energy = self._verify_energy(energy) if energy is not None else self.energy + if include_residuals and custom_energy: + # Residuals are data - model on the experiment grid; mixing them with a model on a + # custom grid would make the DataGroup internally inconsistent. + warnings.warn( + 'Residuals are computed on the experiment energy grid and are omitted ' + 'when a custom energy grid is given.', + UserWarning, + stacklevel=2, + ) + include_residuals = False + data_and_model = { 'Data': self.experiment.binned_data, 'Model': self._create_model_array(energy=energy), @@ -683,8 +697,6 @@ def plot_parameters( plot_kwargs_defaults.update(kwargs) - import plopp as pp - return pp.plot( data_to_plot, **plot_kwargs_defaults, @@ -766,6 +778,18 @@ def _on_convolution_settings_changed(self) -> None: self._owner_index = None self._invalidate_fitter() + def _on_detailed_balance_settings_changed(self) -> None: + """ + Update the detailed balance settings when they change. + + The per-Q analyses hold the settings object they were built with, so replacing it on this + Analysis requires rebuilding the list for the new object to reach every Q index. + """ + super()._on_detailed_balance_settings_changed() + self._analysis_list_is_dirty = True + self._owner_index = None + self._invalidate_fitter() + def _ensure_analysis_list_current(self) -> None: """Rebuild the analysis list if any dependency has changed since it was last built.""" if self._analysis_list_is_dirty and self.Q is not None: @@ -953,17 +977,36 @@ def _fit_all_Q_simultaneously(self) -> FitResults: energy=self.experiment.get_masked_energy(Q_index=analysis1d.Q_index, mask=mask_var) ) - mf = MultiFitter( - fit_objects=self.analysis_list, - fit_functions=self.get_fit_functions(), - ) - - return mf.fit( + # Use the configured fitter rather than a throwaway MultiFitter, so minimizer and + # tolerance settings applied through the ``fitter`` property take effect. + return self.fitter.fit( x=xs, y=ys, weights=ws, ) + def get_all_variables(self) -> list[Parameter]: + """ + Get all variables used in the analysis, across every Q index. + + Overrides the easyscience fallback, which scans every attribute of the object and would + therefore build the MultiFitter and the Sampler as side effects of merely listing variables + (and fail outright on an empty analysis). + + Returns + ------- + list[Parameter] + A list of all variables, including any extra parameters. + """ + variables = self.sample_model.get_all_variables() + + variables.extend(self.instrument_model.get_all_variables()) + + if self._extra_parameters: + variables.extend(self._extra_parameters) + + return variables + def get_fit_functions(self) -> list[callable]: """ Get fit functions for all Q indices, which can be used for simultaneous fitting. @@ -1057,9 +1100,10 @@ def _create_components_dataset( ############# def __repr__(self) -> str: + # The property ensures the list is current, so n_analyses is not reported stale. return ( f'{self.__class__.__name__}(' f'display_name={self.display_name!r}, ' f'unique_name={self.unique_name!r}, ' - f'n_analyses={len(self._analysis_list)})' + f'n_analyses={len(self.analysis_list)})' ) diff --git a/src/easydynamics/analysis/analysis1d.py b/src/easydynamics/analysis/analysis1d.py index 2814e269d..6ac0b7eb3 100644 --- a/src/easydynamics/analysis/analysis1d.py +++ b/src/easydynamics/analysis/analysis1d.py @@ -1,9 +1,11 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import warnings from typing import Any import numpy as np +import plopp as pp import scipp as sc from easyscience.fitting.fitter import Fitter as EasyScienceFitter from easyscience.fitting.minimizers.utils import FitResults @@ -120,6 +122,9 @@ def __init__( self._fit_result = None self._convolver = None self._convolver_is_dirty = True + # The model state_versions the convolver was built against; None until it is built. + # Tracked per Analysis1d so sibling analyses sharing a model each notice a change. + self._convolver_model_versions = None self._fitter = None self._fitter_is_dirty = True self._bayesian = None @@ -360,15 +365,34 @@ def _prepare_for_sampling(self) -> None: The energy grid is fixed for the duration of a fit or a sampling run, so the convolution objects are built once here and reused for every model evaluation. + + Staleness is detected by comparing the models' ``state_version`` against the versions the + convolver was built with. Unlike polling the models' dirty flags, reading a version + consumes nothing, so every Analysis1d sharing a model notices the change — not just the + first one to ask. """ - if ( - self.sample_model.component_collections_is_dirty - or self.instrument_model.resolution_model.component_collections_is_dirty - ): + current = self._model_state_versions() + if None in current or current != self._convolver_model_versions: self._convolver_is_dirty = True self._ensure_convolver_current() + def _model_state_versions(self) -> tuple: + """ + Get the current ``state_version`` of each model the convolver depends on. + + Returns + ------- + tuple + The ``(sample_model, resolution_model)`` state versions. ``None`` entries, for models + that do not expose ``state_version`` yet, never compare equal to a recorded build + version, so the convolver is then conservatively rebuilt. + """ + return ( + getattr(self.sample_model, 'state_version', None), + getattr(self.instrument_model.resolution_model, 'state_version', None), + ) + def as_fit_function( self, _x: np.ndarray | sc.Variable | None = None, @@ -455,8 +479,6 @@ def plot_data_and_model( InteractiveFigure A plot of the data and model. """ - import plopp as pp - data_and_model = self.data_and_model_to_datagroup( energy=energy, add_background=add_background, @@ -467,7 +489,8 @@ def plot_data_and_model( plot_kwargs_defaults = self._build_plot_style_defaults(data_and_model) plot_kwargs_defaults.update(kwargs) - if plot_residuals: + # Residuals may have been omitted (with a warning) for a custom energy grid. + if plot_residuals and 'Residuals' in data_and_model: fig = slicerplot_with_residuals( data_and_model, residuals_key='Residuals', @@ -539,10 +562,22 @@ def data_and_model_to_datagroup( raise ValueError('Q_index must be set to create DataGroup.') energy = self._verify_energy(energy) + custom_energy = energy is not None if energy is None: energy = self._masked_energy + if include_residuals and custom_energy: + # Residuals are data - model on the experiment grid; mixing them with a model on a + # custom grid would make the DataGroup internally inconsistent. + warnings.warn( + 'Residuals are computed on the experiment energy grid and are omitted ' + 'when a custom energy grid is given.', + UserWarning, + stacklevel=2, + ) + include_residuals = False + data_and_model = { 'Data': self.experiment.get_masked_binned_data(Q_index=self.Q_index), 'Model': self._create_model_array(energy=energy), @@ -589,6 +624,7 @@ def rebin(self, dimensions: dict[str, int | sc.Variable]) -> None: def refresh_convolver(self, energy: sc.Variable | None = None) -> None: """Refresh the pre-built Convolution object for the current Q index.""" + self._convolver_model_versions = self._model_state_versions() self._convolver = self._create_convolver(energy=energy) self._convolver_is_dirty = False @@ -660,9 +696,15 @@ def _on_convolution_settings_changed(self) -> None: super()._on_convolution_settings_changed() self._convolver_is_dirty = True + def _on_detailed_balance_settings_changed(self) -> None: + """Mark the convolver as dirty when the detailed balance settings change.""" + super()._on_detailed_balance_settings_changed() + self._convolver_is_dirty = True + def _ensure_convolver_current(self) -> None: """Rebuild the convolver if any dependency has changed since it was last built.""" if self._convolver_is_dirty: + self._convolver_model_versions = self._model_state_versions() self._convolver = self._create_convolver() self._convolver_is_dirty = False diff --git a/src/easydynamics/analysis/analysis_base.py b/src/easydynamics/analysis/analysis_base.py index 48b4abc89..9fe0a531b 100644 --- a/src/easydynamics/analysis/analysis_base.py +++ b/src/easydynamics/analysis/analysis_base.py @@ -375,6 +375,7 @@ def detailed_balance_settings(self, value: DetailedBalanceSettings) -> None: if not isinstance(value, DetailedBalanceSettings): raise TypeError('detailed_balance_settings must be a DetailedBalanceSettings') self._detailed_balance_settings = value + self._on_detailed_balance_settings_changed() @property def extra_parameters(self) -> list[Parameter]: @@ -494,6 +495,11 @@ def _on_convolution_settings_changed(self) -> None: For subclasses that implement convolution, this method can be overridden """ + def _on_detailed_balance_settings_changed(self) -> None: + """ + For subclasses that apply detailed balance, this method can be overridden + """ + def _verify_energy(self, energy: sc.Variable | None) -> sc.Variable | None: """ Verify that the provided energy is the correct type. diff --git a/src/easydynamics/analysis/parameter_analysis.py b/src/easydynamics/analysis/parameter_analysis.py index ffe5fd2cf..1e99f006a 100644 --- a/src/easydynamics/analysis/parameter_analysis.py +++ b/src/easydynamics/analysis/parameter_analysis.py @@ -270,8 +270,8 @@ def fit(self) -> FitResults: The results of the fit """ - xs, ys, ws, _, models = self._build_fit_inputs() - self._invalidate_fitter_if_targets_changed(models) + xs, ys, ws, _, _ = self._build_fit_inputs() + self._invalidate_fitter_if_targets_changed() return self.fitter.fit(x=xs, y=ys, weights=ws) def _build_fit_inputs(self) -> tuple[list, list, list, list, list]: @@ -343,43 +343,42 @@ def _build_fitter(self) -> MultiFitter: A MultiFitter over the per-target models and fit functions. """ _, _, _, funcs, models = self._build_fit_inputs() - self._fitter_targets = self._target_signature(models) + self._fitter_targets = self._target_signature() return MultiFitter(fit_objects=models, fit_functions=funcs) - @staticmethod - def _target_signature(models: list) -> tuple: + def _target_signature(self) -> tuple: """ - Summarize which models the fitter was built for, in target order. + Summarize what the fitter was built for, in target order. - Parameters - ---------- - models : list - The model behind each fit target. + Each entry records the target's model, prediction name, and dataset key. The targets + themselves must be part of the signature, not just the models: swapping which predictions a + binding fits (``binding.targets = ['width'] -> ['area']``) keeps the model list identical + while changing both the frozen fit functions and the data they are fitted against. Returns ------- tuple A comparable signature of the current targets. """ - return tuple(model.unique_name for model in models) + return tuple( + (binding.model.unique_name, target.name, target.dataset_key) + for binding in self.bindings + for target in binding.get_targets() + ) - def _invalidate_fitter_if_targets_changed(self, models: list) -> None: + def _invalidate_fitter_if_targets_changed(self) -> None: """ Rebuild the cached fitter when the bindings no longer resolve to the same targets. A FitBinding can be edited in place -- ``binding.targets = ...`` -- which this object - cannot observe. Doing so changes how many datasets there are, while the cached MultiFitter - still holds the old fit functions, and the fit then dies deep inside the minimizer. Compare - the targets the fitter was built for against the current ones instead. - - Parameters - ---------- - models : list - The model behind each fit target, as currently resolved. + cannot observe. Doing so changes which functions are fitted against which datasets, while + the cached MultiFitter still holds the old fit functions, and the fit then either dies deep + inside the minimizer or silently fits stale functions. Compare the targets the fitter was + built for against the current ones instead. """ if self._fitter is None: return - if self._target_signature(models) != getattr(self, '_fitter_targets', None): + if self._target_signature() != getattr(self, '_fitter_targets', None): self._invalidate_fitter() def _sampling_data(self) -> tuple[list, list, list]: @@ -391,8 +390,8 @@ def _sampling_data(self) -> tuple[list, list, list]: tuple[list, list, list] The ``(x, y, weights)`` triple, one entry per fit target. """ - xs, ys, ws, _, models = self._build_fit_inputs() - self._invalidate_fitter_if_targets_changed(models) + xs, ys, ws, _, _ = self._build_fit_inputs() + self._invalidate_fitter_if_targets_changed() return xs, ys, ws def _chain_parameters(self) -> list[Parameter]: @@ -445,10 +444,6 @@ def plot( if self.parameters is None: raise ValueError('No parameters available to plot.') - full_model_dataset = None - if self.bindings: - full_model_dataset = self.calculate_model_dataset(self.bindings) - # If no names are provided, default to plot all parameters that have bindings. # If no bindings are provided, plot all parameters. if names is None: @@ -462,6 +457,23 @@ def plot( names = self._normalize_names(names) + if not names: + raise ValueError( + 'names must not be an empty list. Pass parameter names to plot, ' + 'or None to plot all parameters with bindings.' + ) + + # Evaluate only the bindings whose targets are actually being plotted. + full_model_dataset = None + if self.bindings: + relevant_bindings = [ + b + for b in self.bindings + if any(target.dataset_key in names for target in b.get_targets()) + ] + if relevant_bindings: + full_model_dataset = self.calculate_model_dataset(relevant_bindings) + # Check that the units of the specified parameters are consistent. units = [self.parameters[name].unit for name in names] first_unit = units[0] @@ -589,12 +601,14 @@ def append_binding(self, binding: FitBinding) -> None: if not isinstance(binding, FitBinding): raise TypeError('binding must be a FitBinding object.') self._bindings.append(binding) + self._invalidate_fitter() def clear_bindings(self) -> None: """ Clear all FitBindings from the list of bindings for the parameter analysis. """ self._bindings.clear() + self._invalidate_fitter() def get_all_variables(self) -> list: """ @@ -638,7 +652,8 @@ def _verify_bindings(self, bindings: FitBinding | list[FitBinding] | None) -> li if isinstance(bindings, FitBinding): return [bindings] if isinstance(bindings, list) and all(isinstance(b, FitBinding) for b in bindings): - return bindings + # Copy so later mutation of the caller's list cannot silently change the bindings. + return list(bindings) raise TypeError('bindings must be a FitBinding, a list of FitBindings, or None.') def _verify_parameters(self, parameters: sc.Dataset | Analysis | None) -> sc.Dataset | None: @@ -799,7 +814,10 @@ def _get_xyweight_from_dataset( q_values = self._parameters[parameter_name].coords['Q'].values if variances is None: - return q_values, values, np.ones_like(values) + # Apply the same finite filtering as the variance path: NaN values arise when a + # parameter is absent for a given Q, and must not leak into a fit. + finite_mask = np.isfinite(values) + return q_values[finite_mask], values[finite_mask], np.ones_like(values[finite_mask]) # NaN variances arise when a parameter is absent for a given Q (parameters_to_dataset # fills np.nan for missing parameters). Filter those rows silently; other non-finite or diff --git a/src/easydynamics/analysis/posterior_sampling.py b/src/easydynamics/analysis/posterior_sampling.py index 04ca3104a..7d2615938 100644 --- a/src/easydynamics/analysis/posterior_sampling.py +++ b/src/easydynamics/analysis/posterior_sampling.py @@ -783,6 +783,8 @@ def plot_trace(self, **kwargs: dict[str, Any]) -> Figure: Figure The matplotlib Figure. """ + # Deliberately imported lazily, to guard against an import cycle between the + # analysis and utils packages. from easydynamics.utils.posterior_plotting import plot_trace results = self._require_results() @@ -809,6 +811,8 @@ def plot_corner(self, **kwargs: dict[str, Any]) -> Figure: Figure The matplotlib Figure. """ + # Deliberately imported lazily, to guard against an import cycle between the + # analysis and utils packages. from easydynamics.utils.posterior_plotting import plot_corner results = self._require_results() @@ -909,6 +913,8 @@ def plot_posterior_predictive( ValueError If n_draws is not a positive integer. """ + # Deliberately imported lazily, to guard against an import cycle between the + # analysis and utils packages. from easydynamics.utils.posterior_plotting import plot_posterior_predictive if not isinstance(n_draws, int) or isinstance(n_draws, bool) or n_draws < 1: @@ -1203,8 +1209,8 @@ def sample( Notes ----- - An ``IndexError`` or ``TypeError`` propagates from the Q_index validation if Q_index is - out of range or not an int. + An ``IndexError`` or ``TypeError`` propagates from the Q_index validation if Q_index is out + of range or not an int. """ if fit_method not in ('independent', 'simultaneous'): raise ValueError("Invalid fit method. Choose 'independent' or 'simultaneous'.") @@ -1400,9 +1406,11 @@ def plot_corner(self, Q_index: int | None = None, **kwargs: dict[str, Any]) -> F Notes ----- - An ``IndexError`` or ``TypeError`` propagates from the Q_index validation if Q_index is - out of range or not an int. + An ``IndexError`` or ``TypeError`` propagates from the Q_index validation if Q_index is out + of range or not an int. """ + # Deliberately imported lazily, to guard against an import cycle between the + # analysis and utils packages. from easydynamics.utils.posterior_plotting import corner_with_slider verify_Q_index(Q_index=Q_index, Q=self._analysis.Q, allow_none=True) @@ -1614,9 +1622,9 @@ def plot_posterior_predictive( ----- A ``NotImplementedError`` propagates when the latest chain is simultaneous: it binds every dataset at once, and no per-Q chain exists for Q_index to pick out. A ``RuntimeError`` - propagates if a slider is asked for outside a notebook or nothing has been sampled yet, - and an ``IndexError`` or ``TypeError`` from the Q_index validation if Q_index is out of - range or not an int. + propagates if a slider is asked for outside a notebook or nothing has been sampled yet, and + an ``IndexError`` or ``TypeError`` from the Q_index validation if Q_index is out of range + or not an int. """ if not isinstance(n_draws, int) or isinstance(n_draws, bool) or n_draws < 1: raise ValueError(f'n_draws must be a positive integer. Got {n_draws}.') diff --git a/src/easydynamics/base_classes/easydynamics_list.py b/src/easydynamics/base_classes/easydynamics_list.py index 74e48a5b9..7c80b2e95 100644 --- a/src/easydynamics/base_classes/easydynamics_list.py +++ b/src/easydynamics/base_classes/easydynamics_list.py @@ -12,7 +12,7 @@ from easydynamics.base_classes.easydynamics_modelbase import EasyDynamicsModelBase from easydynamics.exceptions import AmbiguousNameError -ProtectedType_ = TypeVar('T', bound=EasyDynamicsBase | EasyDynamicsModelBase) +ProtectedType_ = TypeVar('ProtectedType_', bound=EasyDynamicsBase | EasyDynamicsModelBase) class EasyDynamicsList(EasyList[ProtectedType_]): @@ -49,6 +49,10 @@ def __init__( if display_name is None: display_name = unique_name + # Must exist before super().__init__, which appends the initial items through the + # version-bumping mutators below. + self._version = 0 + super().__init__( *args, protected_types=protected_types, @@ -57,6 +61,31 @@ def __init__( **kwargs, ) + # A freshly constructed list always reports version 0, regardless of how many + # initial items were added during construction. + self._version = 0 + + # ------------------------------------------------------------------ + # Properties + # ------------------------------------------------------------------ + + @property + def version(self) -> int: + """ + Get the mutation version of the list. + + Starts at 0 for a freshly constructed list and is incremented by every mutating operation + (append, insert, extend, remove, pop, clear, sort, item assignment and deletion). Consumers + can record the version and later compare it to detect in-place mutations without callbacks. + Read-only; reading never mutates the list. + + Returns + ------- + int + The current mutation version. + """ + return self._version + # ------------------------------------------------------------------ # List methods # ------------------------------------------------------------------ @@ -87,6 +116,7 @@ def insert(self, index: int, value: ProtectedType_) -> None: return super().insert(index, value) + self._bump_version() def append(self, value: ProtectedType_) -> None: """ @@ -126,14 +156,32 @@ def pop(self, index: int | str = -1) -> ProtectedType_: # Overwritten to update warning if isinstance(index, int): - return self._data.pop(index) + item = self._data.pop(index) + self._bump_version() + return item if isinstance(index, str): for i, item in enumerate(self._data): if self._get_key(item) == index: - return self._data.pop(i) + popped = self._data.pop(i) + self._bump_version() + return popped raise KeyError(f'No item with name "{index}" found') raise TypeError('Index must be an int or str') + def sort(self, key: object = None, reverse: bool = False) -> None: + """ + Sort the list in place according to the given key function. + + Parameters + ---------- + key : object, default=None + Mapping function to sort by. + reverse : bool, default=False + Whether to reverse the sort order. + """ + super().sort(key=key, reverse=reverse) + self._bump_version() + # ------------------------------------------------------------------ # Other methods # ------------------------------------------------------------------ @@ -172,6 +220,30 @@ def get_duplicate_names(self) -> list[str]: # Private methods # ------------------------------------------------------------------ + def _bump_version(self) -> None: + """Record that the list was mutated, so version-based consumers rebuild.""" + self._version += 1 + + def _copy_with_items(self, items: list[ProtectedType_]) -> EasyDynamicsList[ProtectedType_]: + """ + Create a new instance of this list class containing the given items. + + Used by slicing. Subclasses whose constructor signature differs from EasyDynamicsList's + (e.g. ComponentCollection) must override this so slicing returns a working instance of the + same class. + + Parameters + ---------- + items : list[ProtectedType_] + The items the new list should contain. + + Returns + ------- + EasyDynamicsList[ProtectedType_] + A new list of the same class containing the items. + """ + return self.__class__(items, protected_types=self._protected_types) + def _get_key(self, obj: EasyDynamicsBase | EasyDynamicsModelBase) -> str: """ Get the name of an object. @@ -241,7 +313,7 @@ def __getitem__( if isinstance(idx, int): return self._data[idx] if isinstance(idx, slice): - return self.__class__(self._data[idx], protected_types=self._protected_types) + return self._copy_with_items(self._data[idx]) if isinstance(idx, str): matches = [r for r in self._data if self._get_key(r) == idx] if len(matches) == 1: @@ -251,3 +323,57 @@ def __getitem__( raise KeyError(f'No item with name "{idx}" found') raise TypeError('Index must be an int, slice, or str') + + def __setitem__(self, idx: int | slice, value: ProtectedType_ | list[ProtectedType_]) -> None: + """ + Set an item (or slice of items) in the list. + + Mirrors the duplicate handling of append/insert: assigning an item that is already in the + list (to a different position) warns and is ignored. + + Parameters + ---------- + idx : int | slice + The index or slice to assign to. + value : ProtectedType_ | list[ProtectedType_] + The new item (or items, for a slice) to assign. Items must be instances of one of the + protected types. + + Notes + ----- + A ``TypeError`` propagates from the type validation or the base assignment if idx or value + has an invalid type, and a ``ValueError`` propagates from the base assignment if slice + assignment changes the slice length. + """ + if isinstance(idx, int): + self._validate_type(value) + if value is not self._data[idx] and value in self: + warnings.warn( + ( + f'Item with name "{self._get_key(value)}" already ' + f'in EasyDynamicsList, it will be ignored' + ), + UserWarning, + stacklevel=2, + ) + return + + super().__setitem__(idx, value) + self._bump_version() + + def __delitem__(self, idx: int | slice | str) -> None: + """ + Delete an item by index, slice, or name. + + Parameters + ---------- + idx : int | slice | str + Index, slice, or name of the item to delete. + + Notes + ----- + A ``KeyError`` propagates from the base deletion if idx is a string that does not match any + item, and a ``TypeError`` propagates from it if idx is not an int, slice, or string. + """ + super().__delitem__(idx) + self._bump_version() diff --git a/src/easydynamics/base_classes/name_mixin.py b/src/easydynamics/base_classes/name_mixin.py index 608ce561d..52b3472af 100644 --- a/src/easydynamics/base_classes/name_mixin.py +++ b/src/easydynamics/base_classes/name_mixin.py @@ -29,9 +29,11 @@ def __init__( If name is not a string. """ - super().__init__(*args, **kwargs) + # Validate before delegating to the parent class so an invalid name fails fast, + # before the parent registers the object in the global map. if not isinstance(name, str): raise TypeError('Name must be a string.') + super().__init__(*args, **kwargs) self._name = name @property diff --git a/src/easydynamics/convolution/analytical_convolution.py b/src/easydynamics/convolution/analytical_convolution.py index efcdf6f52..1f7e4e7ab 100644 --- a/src/easydynamics/convolution/analytical_convolution.py +++ b/src/easydynamics/convolution/analytical_convolution.py @@ -11,7 +11,6 @@ from easydynamics.sample_model import Gaussian from easydynamics.sample_model import Lorentzian from easydynamics.sample_model import Voigt -from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.sample_model.components.model_component import ModelComponent @@ -20,12 +19,13 @@ class AnalyticalConvolution(ConvolutionBase): Analytical convolution of a ModelComponent or ComponentCollection with a ResolutionModel. Possible analytical convolutions are any combination of delta functions, Gaussians, Lorentzians - and Voigt profiles. + and Voigt profiles. Dispatch is subclass-tolerant: a subclass of e.g. Lorentzian is convolved + with the Lorentzian rules. """ - # Mapping of supported component type pairs to convolution methods. + # Mapping of supported canonical component-type-name pairs to convolution methods. # Delta functions are handled separately. - _CONVOLUTIONS: ClassVar[dict[str, object]] = { + _CONVOLUTIONS: ClassVar[dict[tuple[str, str], str]] = { ('Gaussian', 'Gaussian'): '_convolute_gaussian_gaussian', ('Gaussian', 'Lorentzian'): '_convolute_gaussian_lorentzian', ('Gaussian', 'Voigt'): '_convolute_gaussian_voigt', @@ -34,6 +34,65 @@ class AnalyticalConvolution(ConvolutionBase): ('Voigt', 'Voigt'): '_convolute_voigt_voigt', } + # The analytical base types used to resolve a component (or a subclass of one of them) + # to its canonical dispatch name. + _ANALYTICAL_TYPES: ClassVar[tuple[type[ModelComponent], ...]] = (Gaussian, Lorentzian, Voigt) + + def __init__(self, *args: object, **kwargs: object) -> None: + """ + Initialize the AnalyticalConvolution. + + Accepts the same arguments as ConvolutionBase, but requires sample_components and + resolution_components to be provided. + + Parameters + ---------- + *args : object + Positional arguments passed to ConvolutionBase. + **kwargs : object + Keyword arguments passed to ConvolutionBase. + + Raises + ------ + TypeError + If sample_components or resolution_components is None. + """ + super().__init__(*args, **kwargs) + # ConvolutionBase tolerates None collections, but an analytical convolver cannot + # convolve without both models — fail early with a clear error. + if self._sample_components is None: + raise TypeError( + 'sample_components must be a ComponentCollection or ModelComponent, not None.' + ) + if self._resolution_components is None: + raise TypeError( + 'resolution_components must be a ComponentCollection or ModelComponent, not None.' + ) + + @classmethod + def _canonical_type_name(cls, component: ModelComponent) -> str: + """ + Resolve a component to the canonical analytical type name used for dispatch. + + A subclass of one of the analytical types (Gaussian, Lorentzian, Voigt) resolves to its + base type's name, so subclasses are convolved with the base type's rules. + + Parameters + ---------- + component : ModelComponent + The component to resolve. + + Returns + ------- + str + The canonical type name, or the component's own class name if it is not an analytical + type. + """ + for analytical_type in cls._ANALYTICAL_TYPES: + if isinstance(component, analytical_type): + return analytical_type.__name__ + return type(component).__name__ + def convolution( self, ) -> np.ndarray: @@ -90,8 +149,8 @@ def _convolute_analytic_pair( The convolution of two Voigt profiles results in another Voigt profile, with the Gaussian widths summed in quadrature and the Lorentzian widths summed. - The convolution of a delta function with any component or ComponentCollection results in - the same component or ComponentCollection shifted by the delta center. + The convolution of a delta function with any component results in the same component + shifted by the delta center. All areas are multiplied in the convolution. @@ -127,15 +186,15 @@ def _convolute_analytic_pair( resolution_component, ) - pair = (type(sample_component).__name__, type(resolution_component).__name__) + sample_name = self._canonical_type_name(sample_component) + resolution_name = self._canonical_type_name(resolution_component) + + pair = (sample_name, resolution_name) swapped = False if pair not in self._CONVOLUTIONS: # Try reversing the pair - pair = ( - type(resolution_component).__name__, - type(sample_component).__name__, - ) + pair = (resolution_name, sample_name) swapped = True func_name = self._CONVOLUTIONS.get(pair) @@ -154,26 +213,25 @@ def _convolute_analytic_pair( def _convolute_delta_any( self, sample_component: DeltaFunction, - resolution_components: ComponentCollection | ModelComponent, + resolution_component: ModelComponent, ) -> np.ndarray: """ - Convolution of delta function with any ModelComponent or ComponentCollection results in the - same component or ComponentCollection shifted by the delta center. The areas are - multiplied. + Convolution of a delta function with a resolution component results in the same component + shifted by the delta center. The areas are multiplied. Parameters ---------- sample_component : DeltaFunction - The sample component to be convolved. - resolution_components : ComponentCollection | ModelComponent - The resolution model to convolve with. + The sample delta function to be convolved. + resolution_component : ModelComponent + The resolution component to convolve with. Returns ------- np.ndarray The evaluated convolution values at self.energy. """ - return sample_component.area.value * resolution_components.evaluate( + return sample_component.area.value * resolution_component.evaluate( self.energy_with_offset.values - sample_component.center.value ) diff --git a/src/easydynamics/convolution/convolution.py b/src/easydynamics/convolution/convolution.py index 578df564f..7bf29c64e 100644 --- a/src/easydynamics/convolution/convolution.py +++ b/src/easydynamics/convolution/convolution.py @@ -79,13 +79,14 @@ class Convolution(NumericalConvolutionBase): # needs to be rebuilt. # Note: the public 'energy' property setter always writes to '_energy', so '_energy' alone # is sufficient — listing 'energy' separately would cause a double invalidation. + # In-place mutations of the collections, settings-flag changes, and energy_offset + # rebinds are detected separately via the plan-state snapshot and the settings' plan + # versions (see NumericalConvolutionBase._convolution_plan_is_current). _invalidate_plan_on_change: ClassVar[set[str]] = { '_energy', '_sample_components', '_resolution_components', '_temperature', - '_energy_unit', - '_normalize_detailed_balance', '_detailed_balance_settings', } @@ -223,7 +224,7 @@ def _check_if_pair_is_analytic( Raises ------ - TypeError + ValueError If the resolution component is a DeltaFunction. Returns @@ -233,8 +234,8 @@ def _check_if_pair_is_analytic( """ if isinstance(resolution_component, DeltaFunction): - raise TypeError( - 'resolution components contains delta functions. This is not supported.' + raise ValueError( + 'resolution_components contains delta functions. This is not supported.' ) analytical_types = (Gaussian, Lorentzian, Voigt) @@ -243,11 +244,49 @@ def _check_if_pair_is_analytic( and isinstance(resolution_component, analytical_types) ) + def _prune_plan_object(self, obj: object) -> None: + """ + Remove a plan-internal object from the easyscience global map. + + The plan collections and sub-convolvers are private, per-plan objects recreated on every + rebuild; pruning the previous generation keeps the global map from growing with every + rebuild. + + Parameters + ---------- + obj : object + The object to prune, or None for a no-op. + """ + if obj is not None: + self._global_object.map.prune(obj.unique_name) + def _build_convolution_plan(self) -> None: """ Separate sample model components into analytical pairs, delta functions, and the rest. + + Raises + ------ + ValueError + If the resolution collection is empty or contains a DeltaFunction. """ + if self._resolution_components.is_empty: + raise ValueError( + 'resolution_components is empty. Convolution with an empty resolution ' + 'model is not defined; add at least one resolution component.' + ) + self._validate_no_delta_in_resolution(self._resolution_components) + + # Previous plan collections are recreated below; remove them from the global map so + # rebuilds do not leak registry entries. + self._prune_plan_object(getattr(self, '_analytical_sample_components', None)) + self._prune_plan_object(getattr(self, '_delta_sample_components', None)) + self._prune_plan_object(getattr(self, '_numerical_sample_components', None)) + + # Keep the (otherwise unused) inherited dense grid in sync with the current energy + # and settings so it can never hold stale state. + self._energy_grid = self._create_energy_grid() + analytical_sample_components = ComponentCollection(x_unit=self.x_unit, y_unit=self.y_unit) delta_sample_components = ComponentCollection(x_unit=self.x_unit, y_unit=self.y_unit) numerical_sample_components = ComponentCollection(x_unit=self.x_unit, y_unit=self.y_unit) @@ -299,6 +338,11 @@ def _set_convolvers(self) -> None: convolution method. """ + # Previous sub-convolvers are recreated below; remove them from the global map so + # rebuilds do not leak registry entries. + self._prune_plan_object(getattr(self, '_analytical_convolver', None)) + self._prune_plan_object(getattr(self, '_numerical_convolver', None)) + if self._analytical_sample_components: self._analytical_convolver = AnalyticalConvolution( energy=self.energy, @@ -337,15 +381,22 @@ def convert_y_unit(self, unit: str) -> None: The new y-axis unit. """ super().convert_y_unit(unit) - # The sub-convolvers share this convolver's component objects, which were already - # converted by super(); only their y-unit labels need updating. + # The sub-convolvers and plan collections share this convolver's component objects, + # which were already converted by super(); only their y-unit labels need updating. if getattr(self, '_analytical_convolver', None) is not None: self._analytical_convolver._relabel_y_unit(self.y_unit) # ruff: ignore[private-member-access] if getattr(self, '_numerical_convolver', None) is not None: self._numerical_convolver._relabel_y_unit(self.y_unit) # ruff: ignore[private-member-access] + for collection in ( + getattr(self, '_analytical_sample_components', None), + getattr(self, '_delta_sample_components', None), + getattr(self, '_numerical_sample_components', None), + ): + if collection is not None: + collection._y_unit = self.y_unit # ruff: ignore[private-member-access] # Update some setters so the internal sample models are updated - def __setattr__(self, name: str, value: any) -> None: + def __setattr__(self, name: str, value: object) -> None: """ Custom setattr to invalidate convolution plan on relevant attribute changes, and build a new plan. @@ -357,7 +408,7 @@ def __setattr__(self, name: str, value: any) -> None: ---------- name : str The name of the attribute to set. - value : any + value : object The value to set the attribute to. """ super().__setattr__(name, value) diff --git a/src/easydynamics/convolution/convolution_base.py b/src/easydynamics/convolution/convolution_base.py index 21eec5e10..9ba0cc96c 100644 --- a/src/easydynamics/convolution/convolution_base.py +++ b/src/easydynamics/convolution/convolution_base.py @@ -9,6 +9,7 @@ from easydynamics.base_classes import EasyDynamicsModelBase from easydynamics.sample_model.component_collection import ComponentCollection +from easydynamics.sample_model.components.delta_function import DeltaFunction from easydynamics.sample_model.components.model_component import ModelComponent from easydynamics.utils.utils import Numeric from easydynamics.utils.utils import convert_parameter_unit @@ -59,10 +60,16 @@ def __init__( Raises ------ TypeError - If energy is not a numpy ndarray or a scipp Variable or if energy_unit is not a string - or scipp unit, or if energy_offset is not a number or a Parameter, or if - sample_components is not a ComponentCollection or ModelComponent, or if - resolution_components is not a ComponentCollection or ModelComponent. + If energy is not a numpy ndarray or a scipp Variable or if x_unit is not a string or + scipp unit, or if energy_offset is not a number or a Parameter, or if sample_components + is not a ComponentCollection or ModelComponent, or if resolution_components is not a + ComponentCollection or ModelComponent. + + Notes + ----- + A ``ValueError`` propagates from the validation helpers if resolution_components contains a + DeltaFunction, or if the x_unit of the sample or resolution components does not match the + convolver's x_unit. """ super().__init__( @@ -118,8 +125,72 @@ def __init__( x_unit=resolution_components.x_unit, y_unit=resolution_components.y_unit, ) + self._validate_no_delta_in_resolution(resolution_components) self._resolution_components = resolution_components + self._validate_component_x_units() + + # ------------------------------------------------------------------ + # Validation helpers + # ------------------------------------------------------------------ + + @staticmethod + def _validate_no_delta_in_resolution( + resolution_components: ComponentCollection | None, + ) -> None: + """ + Validate that the resolution collection contains no DeltaFunction components. + + Convolving with a delta function in the resolution is not supported on any path + (analytical, numerical, or delta), so the invariant is enforced when the resolution is + bound to the convolver. + + Parameters + ---------- + resolution_components : ComponentCollection | None + The resolution collection to validate. None is skipped. + + Raises + ------ + ValueError + If resolution_components contains a DeltaFunction. + """ + if resolution_components is None: + return + if any(isinstance(component, DeltaFunction) for component in resolution_components): + raise ValueError( + 'resolution_components contains delta functions. This is not supported.' + ) + + def _validate_component_x_units(self) -> None: + """ + Validate that sample and resolution collections use the convolver's x_unit. + + Components in a different (even compatible) x_unit would be evaluated with raw numbers in + the wrong unit, silently producing wrong results. + + Raises + ------ + ValueError + If a collection's x_unit differs from the convolver's x_unit. + """ + if self._x_unit is None: + return + for label, collection in ( + ('sample_components', self._sample_components), + ('resolution_components', self._resolution_components), + ): + if collection is None or collection.x_unit is None: + continue + # Compare as sc.Unit so unit aliases (e.g. 'ueV' vs 'micro-eV') are not false + # mismatches. + if sc.Unit(str(collection.x_unit)) != sc.Unit(str(self._x_unit)): + raise ValueError( + f'{label} has x_unit {str(collection.x_unit)!r}, which does not match the ' + f'convolver x_unit {str(self._x_unit)!r}. Convert the components with ' + f'convert_x_unit before constructing the convolver.' + ) + @property def energy_offset(self) -> Parameter: """ @@ -192,12 +263,15 @@ def energy(self, energy: np.ndarray | sc.Variable) -> None: Parameters ---------- energy : np.ndarray | sc.Variable - 1D array of energy values where the convolution is evaluated. + 1D array of energy values where the convolution is evaluated. A scipp Variable must + carry the convolver's x_unit; the x_unit itself can only be changed via convert_x_unit. Raises ------ TypeError If energy is not a numpy ndarray or a scipp Variable. + ValueError + If energy is a scipp Variable whose unit differs from the convolver's x_unit. """ if isinstance(energy, Numeric): @@ -210,8 +284,15 @@ def energy(self, energy: np.ndarray | sc.Variable) -> None: self._energy = energy_to_scipp(energy, self._energy.unit) if isinstance(energy, sc.Variable): + # Compare as sc.Unit so unit aliases (e.g. 'ueV' vs 'micro-eV') are not false + # mismatches. + if self._x_unit is not None and energy.unit != sc.Unit(str(self._x_unit)): + raise ValueError( + f'energy has unit {str(energy.unit)!r}, which does not match the convolver ' + f'x_unit {str(self._x_unit)!r}. Use convert_x_unit to change the unit, or ' + f'provide energy in {str(self._x_unit)!r}.' + ) self._energy = energy - self._x_unit = energy.unit def convert_x_unit(self, unit: str | sc.Unit) -> None: """ @@ -237,7 +318,9 @@ def convert_x_unit(self, unit: str | sc.Unit) -> None: old_offset_unit = str(self.energy_offset.unit) def _convert_energy(target_unit: str | sc.Unit) -> None: - self.energy = sc.to_unit(self.energy, target_unit) + # Assign the backing field directly: the public setter rejects unit changes + # (convert_x_unit is the one supported route for those). + self._energy = sc.to_unit(self._energy, target_unit) conversions = [ (_convert_energy, unit, old_x_unit), @@ -249,7 +332,8 @@ def _convert_energy(target_unit: str | sc.Unit) -> None: conversions.append((self.resolution_components.convert_x_unit, unit, old_x_unit)) convert_units_with_rollback(conversions) - self._x_unit = unit + # Keep the str contract for x_unit even when an sc.Unit was passed. + self._x_unit = str(unit) if isinstance(unit, sc.Unit) else unit def convert_y_unit(self, unit: str | sc.Unit) -> None: """ @@ -362,6 +446,11 @@ def resolution_components( ------ TypeError If resolution_components is not a ComponentCollection or ModelComponent. + + Notes + ----- + A ``ValueError`` propagates from the validation helper if resolution_components contains a + DeltaFunction. """ if not isinstance(resolution_components, (ComponentCollection, ModelComponent)): raise TypeError( @@ -374,4 +463,5 @@ def resolution_components( x_unit=resolution_components.x_unit, y_unit=resolution_components.y_unit, ) + self._validate_no_delta_in_resolution(resolution_components) self._resolution_components = resolution_components diff --git a/src/easydynamics/convolution/numerical_convolution_base.py b/src/easydynamics/convolution/numerical_convolution_base.py index e5cadda98..77447b1dc 100644 --- a/src/easydynamics/convolution/numerical_convolution_base.py +++ b/src/easydynamics/convolution/numerical_convolution_base.py @@ -81,8 +81,8 @@ def __init__( Raises ------ TypeError - If temperature is not None, a number, or a Parameter, or if temperature_unit is not a - string or sc.Unit. + If sample_components or resolution_components is None, or if temperature is not None, a + number, or a Parameter, or if temperature_unit is not a string or sc.Unit. """ super().__init__( energy=energy, @@ -95,6 +95,17 @@ def __init__( unique_name=unique_name, ) + # ConvolutionBase tolerates None collections, but numerical convolvers cannot + # convolve without both models — fail early with a clear error. + if self._sample_components is None: + raise TypeError( + 'sample_components must be a ComponentCollection or ModelComponent, not None.' + ) + if self._resolution_components is None: + raise TypeError( + 'resolution_components must be a ComponentCollection or ModelComponent, not None.' + ) + if temperature is not None and not isinstance(temperature, (Numeric, Parameter)): raise TypeError('Temperature must be None, a number or a Parameter.') @@ -126,10 +137,13 @@ def _convolution_plan_is_current(self) -> bool: """ Check whether this convolver's plan is up to date. - Plan validity is tracked per convolver so several convolvers can share one - ConvolutionSettings object: each convolver stores the settings' plan version it last - rebuilt against (None after a convolver-local invalidation such as a new energy grid), and - the settings bump their version whenever an accuracy knob changes. + Plan validity is tracked per convolver so several convolvers can share one settings object: + each convolver stores the plan versions of its ConvolutionSettings and + DetailedBalanceSettings it last rebuilt against (None after a convolver-local invalidation + such as a new energy grid), and the settings bump their versions whenever a knob changes. + In addition, a snapshot of the component collections' mutation versions and the + energy_offset binding is compared, so in-place mutations of a live collection (e.g. + append_component) or rebinding the offset to a new Parameter also invalidate the plan. Returns ------- @@ -139,11 +153,40 @@ def _convolution_plan_is_current(self) -> bool: seen_version = getattr(self, '_plan_seen_version', None) if seen_version is None: return False - return self.convolution_settings._plan_valid_for(seen_version) # ruff: ignore[private-member-access] + if not self.convolution_settings._plan_valid_for(seen_version): # ruff: ignore[private-member-access] + return False + seen_db_version = getattr(self, '_plan_seen_db_version', None) + if not self.detailed_balance_settings._plan_valid_for(seen_db_version): # ruff: ignore[private-member-access] + return False + return getattr(self, '_plan_seen_state', None) == self._plan_state_snapshot() def _mark_convolution_plan_current(self) -> None: """Record that this convolver's plan matches its current state and settings.""" self._plan_seen_version = self.convolution_settings._plan_version # ruff: ignore[private-member-access] + self._plan_seen_db_version = self.detailed_balance_settings._plan_version # ruff: ignore[private-member-access] + self._plan_seen_state = self._plan_state_snapshot() + + def _plan_state_snapshot(self) -> tuple: + """ + Snapshot the mutable state the convolution plan was built from. + + Captures the identity and mutation version of the sample and resolution collections (so + both rebinding and in-place mutation are detected) and the identity of the energy_offset + Parameter (so rebinding to a new Parameter invalidates the plan while numeric assignment + mutating the shared Parameter does not). + + Returns + ------- + tuple + A comparable snapshot of the plan-relevant state. + """ + return ( + id(self._sample_components), + self._sample_components.version, + id(self._resolution_components), + self._resolution_components.version, + id(self._energy_offset), + ) @property def convolution_settings(self) -> ConvolutionSettings: @@ -196,6 +239,22 @@ def energy(self, energy: np.ndarray) -> None: ConvolutionBase.energy.fset(self, energy) self._plan_seen_version = None + def convert_x_unit(self, unit: str | sc.Unit) -> None: + """ + Convert the energy axis, energy_offset, and all components to the specified unit, and + invalidate this convolver's plan. + + The dense grid is rebuilt lazily on the next convolution. Other convolvers sharing the same + ConvolutionSettings are unaffected. + + Parameters + ---------- + unit : str | sc.Unit + The unit of the energy. + """ + super().convert_x_unit(unit) + self._plan_seen_version = None + @property def upsample_factor(self) -> Numeric | None: """ @@ -222,7 +281,7 @@ def upsample_factor(self, factor: Numeric | None) -> None: self.convolution_settings.upsample_factor = factor @property - def extension_factor(self) -> float: + def extension_factor(self) -> float | None: """ Get the extension factor. @@ -231,23 +290,24 @@ def extension_factor(self) -> float: Returns ------- - float - The extension factor. + float | None + The extension factor, or None if unset (only valid while upsample_factor is None). """ return self.convolution_settings.extension_factor @extension_factor.setter - def extension_factor(self, factor: Numeric) -> None: + def extension_factor(self, factor: Numeric | None) -> None: """ Set the extension factor. The extension factor determines how much the energy range is extended on both sides before - convolution. 0.2 means extending by 20% of the original energy span on each side. + convolution. 0.2 means extending by 20% of the original energy span on each side. None is + accepted but requires upsample_factor to be None as well before the next convolution. Parameters ---------- - factor : Numeric + factor : Numeric | None The new extension factor. """ self.convolution_settings.extension_factor = factor @@ -331,6 +391,9 @@ def detailed_balance_settings(self, value: DetailedBalanceSettings) -> None: if not isinstance(value, DetailedBalanceSettings): raise TypeError('detailed_balance_settings must be a DetailedBalanceSettings') self._detailed_balance_settings = value + # Convolver-local invalidation: other convolvers sharing the new settings object are + # unaffected. + self._plan_seen_version = None def _create_energy_grid( self, @@ -352,6 +415,11 @@ def _create_energy_grid( EnergyGrid The dense grid created by upsampling and extending energy. """ + # Validate up front so both the upsampled and the non-upsampled path raise the same + # clear error (a single point has no spacing, so no grid can be built from it). + if len(self.energy.values) < 2: + raise ValueError('Energy array must have at least two points.') + if self.upsample_factor is None: # Check if the array is uniformly spaced. energy_diff = np.diff(self.energy.values) @@ -435,13 +503,22 @@ def _check_width_thresholds( # Handle ComponentCollection or ModelComponent components = model if isinstance(model, ComponentCollection) else [model] + # Cover plain-width components as well as Voigt-style components with separate + # gaussian_width/lorentzian_width parameters. + width_attribute_names = ('width', 'gaussian_width', 'lorentzian_width') + for comp in components: - if hasattr(comp, 'width'): - if comp.width.value > LARGE_WIDTH_THRESHOLD * self._energy_grid.energy_span_dense: + for attribute_name in width_attribute_names: + width_param = getattr(comp, attribute_name, None) + if width_param is None: + continue + width_label = attribute_name.replace('_', ' ') + if width_param.value > LARGE_WIDTH_THRESHOLD * self._energy_grid.energy_span_dense: warnings.warn( ( - f"The width of the {model_name} component '{comp.unique_name}' " - f'({comp.width.value}) is large compared to the span of the input ' + f'The {width_label} of the {model_name} component ' + f"'{comp.unique_name}' " + f'({width_param.value}) is large compared to the span of the input ' f'array ({self._energy_grid.energy_span_dense}). ' f'This may lead to inaccuracies in the convolution. ' f'Increase extension_factor to improve accuracy.' @@ -449,11 +526,12 @@ def _check_width_thresholds( UserWarning, stacklevel=3, ) - if comp.width.value < SMALL_WIDTH_THRESHOLD * self._energy_grid.energy_dense_step: + if width_param.value < SMALL_WIDTH_THRESHOLD * self._energy_grid.energy_dense_step: warnings.warn( ( - f"The width of the {model_name} component '{comp.unique_name}' " - f'({comp.width.value}) is small compared to the spacing of the input ' + f'The {width_label} of the {model_name} component ' + f"'{comp.unique_name}' " + f'({width_param.value}) is small compared to the spacing of the input ' f'array ({self._energy_grid.energy_dense_step}). ' f'This may lead to inaccuracies in the convolution. ' f'Increase upsample_factor to improve accuracy.' diff --git a/src/easydynamics/exceptions.py b/src/easydynamics/exceptions.py index e21f30a29..2298cd8e7 100644 --- a/src/easydynamics/exceptions.py +++ b/src/easydynamics/exceptions.py @@ -3,7 +3,23 @@ class AmbiguousNameError(Exception): - def __init__(self, name: str, matches: list[str]) -> None: + """Raised when a name lookup matches more than one element.""" + + def __init__(self, name: str, matches: list[object]) -> None: + """ + Initialize the AmbiguousNameError. + + Parameters + ---------- + name : str + The ambiguous name that was looked up. + matches : list[object] + The elements whose name matched. The elements' unique names are used in the message so + the matches can be told apart. + """ self.name = name self.matches = matches - super().__init__(f"Ambiguous name '{name}' matches {len(matches)} elements: {matches}") + match_names = [ + match.unique_name if hasattr(match, 'unique_name') else str(match) for match in matches + ] + super().__init__(f"Ambiguous name '{name}' matches {len(matches)} elements: {match_names}") diff --git a/src/easydynamics/experiment/experiment.py b/src/easydynamics/experiment/experiment.py index 5be71f109..058df9abe 100644 --- a/src/easydynamics/experiment/experiment.py +++ b/src/easydynamics/experiment/experiment.py @@ -420,9 +420,6 @@ def rebin(self, dimensions: dict[str, int | sc.Variable]) -> None: ) if isinstance(value, float) and value.is_integer(): # I allow eg. 2.0 as well as 2 value = int(value) - # This line can be removed when scipp resize support - # resizing with coordinates - dimensions[dim] = value if not (isinstance(value, (int, sc.Variable))): raise TypeError( f'Dimension values must be integers or sc.Variable. ' diff --git a/src/easydynamics/sample_model/component_collection.py b/src/easydynamics/sample_model/component_collection.py index adf3bc148..18250080f 100644 --- a/src/easydynamics/sample_model/component_collection.py +++ b/src/easydynamics/sample_model/component_collection.py @@ -12,6 +12,7 @@ from easydynamics.base_classes.easydynamics_list import EasyDynamicsList from easydynamics.base_classes.easydynamics_modelbase import EasyDynamicsModelBase +from easydynamics.exceptions import AmbiguousNameError from easydynamics.sample_model.components.model_component import ModelComponent from easydynamics.utils.fit_target import FitTarget from easydynamics.utils.utils import convert_units_with_rollback @@ -69,7 +70,10 @@ def __init__( Parameters ---------- components : ModelComponent | list[ModelComponent] | None, default=None - Initial model components to add to the ComponentCollection. + Initial model components to add to the ComponentCollection. Components are stored by + reference (not copied), so their Parameters stay shared with the objects passed in; + mutating a component mutates it everywhere it is used. Pass a copy if independent + parameters are needed. x_unit : str | sc.Unit, default='meV' Unit of the x-axis (energy, Q, etc.). y_unit : str | sc.Unit, default='dimensionless' @@ -224,11 +228,15 @@ def append_component(self, component: ModelComponent | ComponentCollection) -> N Append a model component or the components from another ComponentCollection to this ComponentCollection. + Components are appended by reference (not copied): their Parameters stay shared with the + passed-in objects, so a fit through one collection updates the same Parameters seen by any + other holder of the component. Pass a copy if independent parameters are needed. + Parameters ---------- component : ModelComponent | ComponentCollection The component to append. If a ComponentCollection is provided, all of its components - will be appended. + will be appended (also by reference). """ if isinstance(component, ComponentCollection): self.extend(component) @@ -281,8 +289,8 @@ def normalize_area(self) -> None: Raises ------ ValueError - If there are no components in the model or if the total area is zero or not finite, - which would prevent normalization. + If there are no components in the model, if any component area is negative, or if the + total area is zero, negative or not finite, which would prevent normalization. """ if not self: raise ValueError('No components in the model to normalize.') @@ -307,12 +315,19 @@ def normalize_area(self) -> None: # units normalize correctly. Dividing each value by the total expressed in the # reference unit makes the areas sum to 1 in that unit. reference_unit = str(area_params[0].unit) - total_area_value = sum( - convert_value_unit(p.value, p.unit, reference_unit) for p in area_params - ) + area_values = [convert_value_unit(p.value, p.unit, reference_unit) for p in area_params] + + negative = [p.name for p, value in zip(area_params, area_values, strict=True) if value < 0] + if negative: + raise ValueError( + f'Negative area(s) found for {negative}; cannot normalize. ' + 'Areas must be non-negative for normalization to be meaningful.' + ) + + total_area_value = sum(area_values) - if total_area_value == 0: - raise ValueError('Total area is zero; cannot normalize.') + if total_area_value <= 0: + raise ValueError('Total area is not positive; cannot normalize.') if not np.isfinite(total_area_value): raise ValueError('Total area is not finite; cannot normalize.') @@ -350,18 +365,27 @@ def evaluate( output : str, default='numpy' 'numpy' returns np.ndarray; 'scipp' returns sc.Variable with y_unit. + Raises + ------ + ValueError + If output is not 'numpy' or 'scipp'. + Returns ------- np.ndarray | sc.Variable Evaluated model values. """ if not self: + # Mirror the validation and 1D output shape of the non-empty path. + if output not in ('numpy', 'scipp'): + raise ValueError(f"output must be 'numpy' or 'scipp', got {output!r}") if isinstance(x, (sc.Variable, sc.DataArray)): - values = np.zeros_like(x.values, dtype=float) dim = x.dims[0] if x.dims else 'x' + raw = x.values if x.dims else x.value else: - values = np.zeros_like(x, dtype=float) dim = 'x' + raw = x + values = np.zeros_like(np.atleast_1d(np.asarray(raw, dtype=float)), dtype=float) if output == 'scipp': return sc.array(dims=[dim], values=values, unit=self.y_unit) return values @@ -396,6 +420,8 @@ def evaluate_component( If name is not a string. KeyError If no component with the given name exists in the collection. + AmbiguousNameError + If more than one component with the given name exists in the collection. Returns ------- @@ -409,6 +435,8 @@ def evaluate_component( matches = [comp for comp in self if comp.name == name] if not matches: raise KeyError(f"No component named '{name}' exists.") + if len(matches) > 1: + raise AmbiguousNameError(name, matches) return matches[0].evaluate(x, output=output) def fix_all_parameters(self) -> None: @@ -425,6 +453,30 @@ def free_all_parameters(self) -> None: # Private methods # ------------------------------------------------------------------ + def _copy_with_items(self, items: list[ModelComponent]) -> ComponentCollection: + """ + Create a new collection of this class containing the given components. + + Used by slicing. Overridden because ComponentCollection's constructor signature differs + from EasyDynamicsList's. The new collection carries this collection's units and references + the same component objects (no copies). + + Parameters + ---------- + items : list[ModelComponent] + The components the new collection should contain. + + Returns + ------- + ComponentCollection + A new collection of the same class containing the components. + """ + return self.__class__( + components=list(items), + x_unit=self.x_unit, + y_unit=self.y_unit, + ) + def _warn_if_duplicate_names(self) -> None: """Warn if any two components share the same name.""" names = [c.name for c in self] diff --git a/src/easydynamics/sample_model/components/damped_harmonic_oscillator.py b/src/easydynamics/sample_model/components/damped_harmonic_oscillator.py index 707823f5a..b4dfe9a45 100644 --- a/src/easydynamics/sample_model/components/damped_harmonic_oscillator.py +++ b/src/easydynamics/sample_model/components/damped_harmonic_oscillator.py @@ -130,14 +130,13 @@ def area(self, value: Numeric) -> None: value : Numeric New area value (in current area unit = x_unit * y_unit). - Raises - ------ - TypeError - If *value* is not a numeric type. + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). """ - if not isinstance(value, Numeric): - raise TypeError('area must be a number') - self._area.value = value + self._set_bounded_parameter_value(self._area, value, 'area') @property def center(self) -> Parameter: @@ -203,7 +202,7 @@ def width(self, value: Numeric) -> None: raise TypeError('width must be a number') if float(value) <= 0: raise ValueError('width must be positive') - self._width.value = value + self._set_bounded_parameter_value(self._width, value, 'width') def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: r""" diff --git a/src/easydynamics/sample_model/components/delta_function.py b/src/easydynamics/sample_model/components/delta_function.py index 539f2970f..0cd2cbe9a 100644 --- a/src/easydynamics/sample_model/components/delta_function.py +++ b/src/easydynamics/sample_model/components/delta_function.py @@ -11,7 +11,10 @@ from easydynamics.sample_model.components.model_component import ModelComponent from easydynamics.utils.utils import Numeric -EPSILON = 1e-8 # tolerance for bin-edge comparisons +# Absolute tolerance for deciding whether the center falls inside the x range. It is expressed +# in the unit x is evaluated in (typically meV), so it only serves to absorb floating-point +# noise at the grid edges — it is not a physically meaningful width. +EPSILON = 1e-8 if TYPE_CHECKING: import scipp as sc @@ -121,14 +124,13 @@ def area(self, value: Numeric) -> None: value : Numeric New area value (in current area unit = x_unit * y_unit). - Raises - ------ - TypeError - If *value* is not a numeric type. + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). """ - if not isinstance(value, Numeric): - raise TypeError('area must be a number') - self._area.value = value + self._set_bounded_parameter_value(self._area, value, 'area') @property def center(self) -> Parameter: @@ -184,12 +186,25 @@ def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndar Zero everywhere, with a single non-zero bin nearest the center when center falls within the x range. + Raises + ------ + ValueError + If x_vals contains a single point. A delta function's evaluated height is ``area / + bin_width``, and a single point defines no bin width. + Notes ----- When ``center`` falls within the x range, the bin nearest to ``center`` receives ``area / bin_width`` rather than zero. In convolutions, the DeltaFunction acts as an identity element (handled by the Convolution class). """ + if x_vals.size == 1: + raise ValueError( + 'A DeltaFunction cannot be evaluated at a single x value: its evaluated height ' + 'is area / bin_width, and a single point defines no bin width. Evaluate on a ' + 'grid of at least two x values.' + ) + center = self._resolve_param_value(self._center, eval_unit) area = self._resolve_param_value(self._area, self._eval_area_unit(eval_unit)) @@ -205,14 +220,11 @@ def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndar i = np.argmin(np.abs(x_sorted - center)) # left half-width - if i == 0: - left = x_sorted[1] - x_sorted[0] if x_sorted.size > 1 else 0.5 - else: - left = x_sorted[i] - x_sorted[i - 1] + left = x_sorted[i] - x_sorted[i - 1] if i > 0 else x_sorted[1] - x_sorted[0] # right half-width if i == x_sorted.size - 1: - right = x_sorted[-1] - x_sorted[-2] if x_sorted.size > 1 else 0.5 + right = x_sorted[-1] - x_sorted[-2] else: right = x_sorted[i + 1] - x_sorted[i] diff --git a/src/easydynamics/sample_model/components/expression_component.py b/src/easydynamics/sample_model/components/expression_component.py index 25e261049..6b1927804 100644 --- a/src/easydynamics/sample_model/components/expression_component.py +++ b/src/easydynamics/sample_model/components/expression_component.py @@ -29,9 +29,15 @@ class ExpressionComponent(ModelComponent): Model component defined by a symbolic expression. The expression must contain ``x`` as the independent variable. All other symbols are treated as - free parameters, which can be accessed and set as attributes after construction. Supported - functions include ``exp``, ``sin``, ``cos``, ``sqrt``, ``erf``, and others — see the - ``_ALLOWED_FUNCS`` class variable for the full list. + free parameters, which can be accessed and set as attributes after construction. Symbol names + that collide with an existing attribute of the class (e.g. ``name`` or ``evaluate``) are + rejected at construction. Supported functions include ``exp``, ``sin``, ``cos``, ``sqrt``, + ``erf``, and others — see the ``_ALLOWED_FUNCS`` class variable for the full list. + + .. warning:: + The expression string is parsed with ``sympy.sympify``, which evaluates the string and + can execute arbitrary code. Only pass expression strings from a trusted source — never + feed it unsanitized user input. Examples -------- @@ -66,7 +72,10 @@ class ExpressionComponent(ModelComponent): Parameters are dimensionless by default. Units can be given per parameter at construction, or relabelled later with ``set_unit`` (the numeric value is kept as-is). When units are in use, the unit of the evaluated expression is derived from the parameter units and x_unit (see - ``output_unit``), and a warning is issued if it does not match y_unit: + ``output_unit``). A derived unit that differs from y_unit but is convertible to it is handled + automatically: the expression is evaluated in a coherent (SI) scale, so parameter units of + mixed scales combine correctly, and the result is expressed in y_unit. A warning is issued only + when the derived unit is dimensionally incompatible with y_unit: ```python expr = edyn.ExpressionComponent( 'A * exp(-(x - x0)**2 / (2*sigma**2))', @@ -166,15 +175,19 @@ def __init__( expression : str The symbolic expression as a string. Must contain 'x' as the independent variable. The symbols ``hbar`` and ``kb`` are provided automatically as read-only physical constants - (in meV*s and meV/K respectively) unless overridden via *parameters*. + (in meV*s and meV/K respectively) unless overridden via *parameters*. The string is + parsed with ``sympy.sympify``, which can execute arbitrary code — only use expression + strings from a trusted source. Symbol names that collide with an existing attribute of + the class (e.g. ``name``, ``evaluate``) are rejected. parameters : dict[str, Numeric] | None, default=None Dictionary of parameter names and their initial values. Parameters that are not given a unit are dimensionless. parameter_units : dict[str, str | sc.Unit] | None, default=None Optional units per parameter name. Each entry sets the unit of the named parameter without rescaling its value (see :meth:`set_unit`), and takes precedence over the unit - of a Parameter instance given in *parameters*. When units are in use, a warning is - issued if the expression's output unit does not match y_unit. + of a Parameter instance given in *parameters*. When units are in use, an output unit + convertible to y_unit rescales the evaluated values into y_unit; a warning is issued + only if the output unit is incompatible with y_unit. x_unit : str | sc.Unit, default='meV' Unit of the x-axis. y_unit : str | sc.Unit, default='dimensionless' @@ -189,8 +202,9 @@ def __init__( Raises ------ ValueError - If the expression is invalid or does not contain 'x', or if parameter_units names a - parameter that is not in the expression. + If the expression is invalid or does not contain 'x', if a symbol name collides with an + existing attribute of the class, or if parameter_units names a parameter that is not in + the expression. TypeError If any parameter value is not numeric, or if parameter_units is not a dictionary. """ @@ -267,6 +281,16 @@ def __init__( if name in self._RESERVED_NAMES: continue + # A symbol shadowing an existing attribute (e.g. 'name', 'evaluate', 'x_unit') + # would silently diverge: reads resolve to the class attribute (since __getattr__ + # only fires when normal lookup fails) while writes hit the parameter. Reject it. + if hasattr(type(self), name) or name in self.__dict__: + raise ValueError( + f"Symbol '{name}' in the expression collides with an existing attribute " + f'of {type(self).__name__}; it could not be accessed as a parameter. ' + f'Rename the symbol in the expression.' + ) + # Physical constants are provided automatically, unless the user explicitly # supplies a parameter with the same name. if name in self._PHYSICAL_CONSTANTS and name not in parameters: @@ -389,16 +413,32 @@ def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndar f'convert x to {self.x_unit} before evaluating.' ) + # When the derived output unit is convertible to y_unit, evaluate in the coherent SI + # scale: every symbol's value is scaled by its unit's SI multiplier, so mixed-scale + # parameter units combine correctly even inside sums (e.g. 1 + D*x**2*tau with D in m^2/s, + # x in 1/angstrom and tau in ps), and the result is expressed in y_unit. Scale-homogeneous + # expressions give the same numbers either way. + scale_into_si = self._output_converts_to_y_unit() + args = [] for name in self._symbol_names: if name == 'x': - args.append(x_vals) + value = x_vals + unit = self.x_unit elif name in self._constants: - args.append(self._constants[name].value) + value = self._constants[name].value + unit = self._constants[name].unit else: - args.append(self._parameters[name].value) + value = self._parameters[name].value + unit = self._parameters[name].unit + if scale_into_si and unit is not None: + value = value * self._si_multiplier(unit) + args.append(value) - return self._func(*args) + result = self._func(*args) + if scale_into_si: + result = result / self._si_multiplier(self.y_unit or 'dimensionless') + return result def get_all_variables(self) -> list[Parameter]: """ @@ -417,8 +457,9 @@ def set_unit(self, name: str, unit: str | sc.Unit) -> None: This relabels the unit: the numeric value, bounds, and variance are kept as-is. Use ``Parameter.convert_unit`` instead to rescale a value into a compatible unit. Issues a - warning if the resulting output unit of the expression no longer matches y_unit. Raises the - same exceptions as :meth:`_relabel_parameter_unit` on invalid input. + warning if the resulting output unit of the expression is incompatible with y_unit (a + convertible output unit rescales evaluated values into y_unit instead). Raises the same + exceptions as :meth:`_relabel_parameter_unit` on invalid input. Parameters ---------- @@ -661,18 +702,75 @@ def _propagate_unit(self, node: sp.Basic) -> sc.Unit: f'Cannot determine units for expression node {node} of type {type(node).__name__}.' ) + def _units_in_use(self) -> bool: + """ + Whether the expression carries unit information at all. + + Returns + ------- + bool + True when the expression uses physical constants or any parameter has a unit other than + dimensionless. Unit-agnostic expressions (all parameters dimensionless) evaluate + without any unit handling. + """ + return bool(self._constants) or any( + str(parameter.unit) != 'dimensionless' for parameter in self._parameters.values() + ) + + @staticmethod + def _si_multiplier(unit: str | sc.Unit) -> float: + """ + Scale factor from one of *unit* to the coherent SI value of the same dimension. + + Parameters + ---------- + unit : str | sc.Unit + The unit whose scale to extract, e.g. 1e-10 for angstrom. + + Returns + ------- + float + The multiplier relative to the coherent SI base units. + """ + return float(sc.Unit(str(unit)).to_dict().get('multiplier', 1.0)) + + def _output_converts_to_y_unit(self) -> bool: + """ + Whether evaluation should run in a coherent scale and express the result in y_unit. + + Returns + ------- + bool + True when units are in use and the derived output unit differs from y_unit but is + convertible to it. False when units are not in use, the output unit cannot be + determined, the units already agree (no conversion needed), or they are dimensionally + incompatible (construction warned; values are evaluated raw and labelled as-is). + """ + if not self._units_in_use(): + return False + try: + output_unit = sc.Unit(self.output_unit) + except sc.UnitError: + return False + y_unit = sc.Unit(self.y_unit) if self.y_unit is not None else sc.Unit('dimensionless') + if output_unit == y_unit: + return False + try: + sc.to_unit(sc.scalar(1.0, unit=output_unit), y_unit) + except sc.UnitError: + return False + return True + def _warn_if_output_unit_mismatch(self) -> None: """ - Warn if the expression's output unit does not match y_unit. + Warn if the expression's output unit cannot be expressed in y_unit. The check only runs when units are in use, i.e. when the expression uses physical constants or any parameter has a unit other than dimensionless. Unit-agnostic expressions (all - parameters dimensionless) stay silent. + parameters dimensionless) stay silent. An output unit that differs from y_unit but is + convertible to it does not warn: evaluated values are rescaled into y_unit. """ - units_in_use = bool(self._constants) or any( - str(parameter.unit) != 'dimensionless' for parameter in self._parameters.values() - ) - if not units_in_use: + if not self._units_in_use(): return try: @@ -686,11 +784,16 @@ def _warn_if_output_unit_mismatch(self) -> None: return y_unit = sc.Unit(self.y_unit) if self.y_unit is not None else sc.Unit('dimensionless') - if output_unit != y_unit: + if output_unit == y_unit: + return + try: + sc.to_unit(sc.scalar(1.0, unit=output_unit), y_unit) + except sc.UnitError: warnings.warn( f'The expression evaluates to unit {output_unit}, which does not match ' - f'y_unit {y_unit}. The evaluated values are labelled with y_unit; adjust the ' - f'parameter units or y_unit to make them consistent.', + f'y_unit {y_unit} and cannot be converted to it. The evaluated values are ' + f'labelled with y_unit; adjust the parameter units or y_unit to make them ' + f'consistent.', UserWarning, stacklevel=3, ) diff --git a/src/easydynamics/sample_model/components/gaussian.py b/src/easydynamics/sample_model/components/gaussian.py index db9ad082a..65a09a528 100644 --- a/src/easydynamics/sample_model/components/gaussian.py +++ b/src/easydynamics/sample_model/components/gaussian.py @@ -126,14 +126,13 @@ def area(self, value: Numeric) -> None: value : Numeric New area value (in current area unit = x_unit * y_unit). - Raises - ------ - TypeError - If *value* is not a numeric type. + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). """ - if not isinstance(value, Numeric): - raise TypeError('area must be a number') - self._area.value = value + self._set_bounded_parameter_value(self._area, value, 'area') @property def center(self) -> Parameter: @@ -193,13 +192,13 @@ def width(self, value: Numeric) -> None: TypeError If *value* is not a numeric type. ValueError - If *value* is not positive. + If *value* is not positive, or violates the width parameter's bounds. """ if not isinstance(value, Numeric): raise TypeError('width must be a number') if float(value) <= 0: raise ValueError('width must be positive') - self._width.value = value + self._set_bounded_parameter_value(self._width, value, 'width') def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: r""" diff --git a/src/easydynamics/sample_model/components/lorentzian.py b/src/easydynamics/sample_model/components/lorentzian.py index 6aade54f2..5e9ccd74b 100644 --- a/src/easydynamics/sample_model/components/lorentzian.py +++ b/src/easydynamics/sample_model/components/lorentzian.py @@ -22,7 +22,7 @@ class Lorentzian(CreateParametersMixin, ModelComponent): $$ I(x) = \frac{A}{\pi} \frac{\Gamma}{(x - x_0)^2 + \Gamma^2} $$ - where $A$ is the area, $x_0$ is the center, and $\Gamma$ is the hald width at half max (HWHM). + where $A$ is the area, $x_0$ is the center, and $\Gamma$ is the half width at half max (HWHM). area has unit = x_unit * y_unit; center and width have unit = x_unit. If the center is not provided, it will be centered at 0 and fixed, which is typically what you @@ -124,14 +124,13 @@ def area(self, value: Numeric) -> None: value : Numeric New area value (in current area unit = x_unit * y_unit). - Raises - ------ - TypeError - If *value* is not a numeric type. + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). """ - if not isinstance(value, Numeric): - raise TypeError('area must be a number') - self._area.value = value + self._set_bounded_parameter_value(self._area, value, 'area') @property def center(self) -> Parameter: @@ -191,13 +190,13 @@ def width(self, value: Numeric) -> None: TypeError If *value* is not a numeric type. ValueError - If *value* is not positive. + If *value* is not positive, or violates the width parameter's bounds. """ if not isinstance(value, Numeric): raise TypeError('width must be a number') if float(value) <= 0: raise ValueError('width must be positive') - self._width.value = value + self._set_bounded_parameter_value(self._width, value, 'width') def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: r""" diff --git a/src/easydynamics/sample_model/components/mixins.py b/src/easydynamics/sample_model/components/mixins.py index 2fc0b8071..531ab0091 100644 --- a/src/easydynamics/sample_model/components/mixins.py +++ b/src/easydynamics/sample_model/components/mixins.py @@ -21,6 +21,42 @@ class CreateParametersMixin: area_unit = x_unit * y_unit, so when y_unit='dimensionless', area_unit = x_unit. """ + @staticmethod + def _set_bounded_parameter_value(param: Parameter, value: Numeric, label: str) -> None: + """ + Assign a value to a bounded parameter, raising instead of silently clamping. + + easyscience's ``Parameter.value`` setter silently clamps out-of-bounds values to the + nearest bound, which corrupts the parameter (e.g. assigning -1.0 to an area with ``min=0`` + stores 0.0). Component setters route assignments through this helper so a bounds violation + raises a clear error instead. + + Parameters + ---------- + param : Parameter + The parameter to assign to. + value : Numeric + The new value. + label : str + Name of the parameter used in error messages (e.g. ``'area'``, ``'width'``). + + Raises + ------ + TypeError + If *value* is not a numeric type. + ValueError + If *value* violates the parameter's bounds. + """ + if not isinstance(value, Numeric): + raise TypeError(f'{label} must be a number') + value = float(value) + if value < param.min or value > param.max: + raise ValueError( + f'Cannot set {label} to {value}: it violates the parameter bounds ' + f'[{param.min}, {param.max}]. Adjust the bounds first if this value is intended.' + ) + param.value = value + def _create_area_parameter( self, area: Numeric, diff --git a/src/easydynamics/sample_model/components/polynomial.py b/src/easydynamics/sample_model/components/polynomial.py index 4e592cf62..5bcd6f474 100644 --- a/src/easydynamics/sample_model/components/polynomial.py +++ b/src/easydynamics/sample_model/components/polynomial.py @@ -5,6 +5,7 @@ import warnings from collections.abc import Sequence +from contextlib import suppress import numpy as np import scipp as sc @@ -379,12 +380,39 @@ def get_all_variables(self) -> list[DescriptorBase]: """ return list(self._coefficients) + @staticmethod + def _rescale_coefficient(param: Parameter, factor: float) -> None: + """ + Rescale a coefficient's value and bounds by a positive factor without clamping. + + The bounds are temporarily widened to infinity so the value assignment cannot be silently + clamped by the Parameter's min/max (easyscience clamps out-of-bounds values instead of + raising), then the original bounds are rescaled by the same factor. + + Parameters + ---------- + param : Parameter + The coefficient Parameter to rescale. + factor : float + The (strictly positive) rescaling factor. + """ + old_min = param.min + old_max = param.max + param.min = -np.inf + param.max = np.inf + param.value = param.value * factor + param.min = old_min * factor + param.max = old_max * factor + def convert_x_unit(self, new_x_unit: str | sc.Unit) -> None: """ Convert the x-axis unit by rescaling coefficients with power-law factors. Each coefficient ``c_i`` is rescaled by ``(old_scale / new_scale) ** i`` so the evaluated - polynomial output is unchanged after the conversion. + polynomial output is unchanged after the conversion. The coefficient bounds (min/max) are + rescaled by the same factor, so bounded coefficients convert without being clamped. If any + step fails, the already-converted coefficients are rolled back best-effort before the + exception propagates. Parameters ---------- @@ -400,21 +428,37 @@ def convert_x_unit(self, new_x_unit: str | sc.Unit) -> None: if not isinstance(new_x_unit, (str, sc.Unit)): raise UnitError('new_x_unit must be a string or a scipp unit.') - conversion_value_before = self._x_unit_helper.value - self._x_unit_helper = sc.to_unit(self._x_unit_helper, unit=new_x_unit) - conversion_value_after = self._x_unit_helper.value - for i, param in enumerate(self._coefficients): - param.value *= (conversion_value_before / conversion_value_after) ** i - + new_helper = sc.to_unit(self._x_unit_helper, unit=new_x_unit) + scale = self._x_unit_helper.value / new_helper.value + + rescaled: list[tuple[Parameter, float]] = [] + converted = False + try: + for i, param in enumerate(self._coefficients): + factor = scale**i + # Exact comparison on purpose: only a factor of exactly 1.0 (same unit, or the + # constant term's scale**0) is a guaranteed no-op worth skipping. + if factor != 1.0: # ruff: ignore[float-equality-comparison] + self._rescale_coefficient(param, factor) + rescaled.append((param, factor)) + converted = True + finally: + if not converted: + for param, factor in rescaled: + with suppress(Exception): + self._rescale_coefficient(param, 1.0 / factor) + + self._x_unit_helper = new_helper self._x_unit = str(new_x_unit) if isinstance(new_x_unit, sc.Unit) else new_x_unit def convert_y_unit(self, new_y_unit: str | sc.Unit) -> None: """ Rescale all coefficients so the evaluated output remains the same physical value. - All coefficients are multiplied by the conversion factor from ``old_y_unit`` to - ``new_y_unit`` so that ``I(x) [new_y_unit]`` represents the same physical quantity as - ``I(x) [old_y_unit]``. + All coefficients (values and bounds) are multiplied by the conversion factor from + ``old_y_unit`` to ``new_y_unit`` so that ``I(x) [new_y_unit]`` represents the same physical + quantity as ``I(x) [old_y_unit]``. If any step fails, the already-converted coefficients + are rolled back best-effort before the exception propagates. Parameters ---------- @@ -438,8 +482,21 @@ def convert_y_unit(self, new_y_unit: str | sc.Unit) -> None: y_helper_new = sc.to_unit(y_helper, new_y_str) scale = y_helper_new.value / y_helper.value - for param in self._coefficients: - param.value *= scale + rescaled: list[Parameter] = [] + converted = False + try: + for param in self._coefficients: + # Exact comparison on purpose: only a scale of exactly 1.0 (converting to the + # same unit) is a guaranteed no-op worth skipping. + if scale != 1.0: # ruff: ignore[float-equality-comparison] + self._rescale_coefficient(param, scale) + rescaled.append(param) + converted = True + finally: + if not converted: + for param in rescaled: + with suppress(Exception): + self._rescale_coefficient(param, 1.0 / scale) self._y_unit = new_y_str def __repr__(self) -> str: diff --git a/src/easydynamics/sample_model/components/voigt.py b/src/easydynamics/sample_model/components/voigt.py index d0bd070d7..4a196ea7c 100644 --- a/src/easydynamics/sample_model/components/voigt.py +++ b/src/easydynamics/sample_model/components/voigt.py @@ -57,10 +57,10 @@ class Voigt(CreateParametersMixin, ModelComponent): def __init__( self, - area: Numeric | Parameter = 1.0, - center: Numeric | Parameter | None = None, - gaussian_width: Numeric | Parameter = 1.0, - lorentzian_width: Numeric | Parameter = 1.0, + area: Numeric = 1.0, + center: Numeric | None = None, + gaussian_width: Numeric = 1.0, + lorentzian_width: Numeric = 1.0, x_unit: str | sc.Unit = 'meV', y_unit: str | sc.Unit = 'dimensionless', name: str = 'Voigt', @@ -72,13 +72,13 @@ def __init__( Parameters ---------- - area : Numeric | Parameter, default=1.0 + area : Numeric, default=1.0 Integrated area under the Voigt profile. Unit is ``x_unit * y_unit``. - center : Numeric | Parameter | None, default=None + center : Numeric | None, default=None Peak position in x_unit. If None, defaults to 0 and the center parameter is fixed. - gaussian_width : Numeric | Parameter, default=1.0 + gaussian_width : Numeric, default=1.0 Gaussian component standard deviation (sigma) in x_unit. Must be strictly positive. - lorentzian_width : Numeric | Parameter, default=1.0 + lorentzian_width : Numeric, default=1.0 Lorentzian component HWHM (gamma) in x_unit. Must be strictly positive. x_unit : str | sc.Unit, default='meV' Unit of the x-axis. center, gaussian_width, and lorentzian_width are stored in this @@ -139,14 +139,13 @@ def area(self, value: Numeric) -> None: value : Numeric New area value (in current area unit = x_unit * y_unit). - Raises - ------ - TypeError - If *value* is not a numeric type. + Notes + ----- + A ``TypeError`` propagates from the shared value setter if *value* is not a numeric type, + and a ``ValueError`` propagates from it if *value* violates the area parameter's bounds + (e.g. a negative value when the area was created non-negative, giving it ``min=0``). """ - if not isinstance(value, Numeric): - raise TypeError('area must be a number') - self._area.value = value + self._set_bounded_parameter_value(self._area, value, 'area') @property def center(self) -> Parameter: @@ -213,7 +212,7 @@ def gaussian_width(self, value: Numeric) -> None: raise TypeError('gaussian_width must be a number') if float(value) <= 0: raise ValueError('gaussian_width must be positive') - self._gaussian_width.value = value + self._set_bounded_parameter_value(self._gaussian_width, value, 'gaussian_width') @property def lorentzian_width(self) -> Parameter: @@ -247,7 +246,7 @@ def lorentzian_width(self, value: Numeric) -> None: raise TypeError('lorentzian_width must be a number') if float(value) <= 0: raise ValueError('lorentzian_width must be positive') - self._lorentzian_width.value = value + self._set_bounded_parameter_value(self._lorentzian_width, value, 'lorentzian_width') def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndarray: """ diff --git a/src/easydynamics/sample_model/diffusion_model/__init__.py b/src/easydynamics/sample_model/diffusion_model/__init__.py index ceee6588b..778abc129 100644 --- a/src/easydynamics/sample_model/diffusion_model/__init__.py +++ b/src/easydynamics/sample_model/diffusion_model/__init__.py @@ -4,11 +4,13 @@ from easydynamics.sample_model.diffusion_model.brownian_translational_diffusion import ( BrownianTranslationalDiffusion, ) +from easydynamics.sample_model.diffusion_model.delta_lorentz import DeltaLorentz from easydynamics.sample_model.diffusion_model.jump_translational_diffusion import ( JumpTranslationalDiffusion, ) __all__ = [ 'BrownianTranslationalDiffusion', + 'DeltaLorentz', 'JumpTranslationalDiffusion', ] diff --git a/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py b/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py index cfcefae92..c820a37f6 100644 --- a/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py +++ b/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py @@ -237,6 +237,9 @@ def create_component_collections( Create ComponentCollection components for the Brownian translational diffusion model at given Q values. + The created collections are installed on the model (they become the collections returned by + ``get_component_collections``), so the returned list is the live one. + Returns ------- list[ComponentCollection] @@ -295,7 +298,8 @@ def create_component_collections( component_collection_list[i].append_component(lorentzian_component) - return component_collection_list + self._component_collections = component_collection_list + return self._component_collections # ------------------------------------------------------------------ # Private methods diff --git a/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py b/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py index 2e255d424..0424b5f98 100644 --- a/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py +++ b/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py @@ -179,19 +179,10 @@ def __init__( # -------------------------------------------------------------- self._allow_Q_variation = self._create_Q_variation_dict(allow_Q_variation) - self._A_0_list = [] - self._A_1_list = [] - self._lorentzian_width_list = [] - if self.Q is not None: - if self._allow_Q_variation['A_0'] is True: - self._A_0_list, self._A_1_list = self._create_A0_A1_parameter_lists(self.A_0) - - if self._allow_Q_variation['lorentzian_width'] is True: - self._lorentzian_width_list = self._create_lorentzian_width_parameter_list( - self.lorentzian_width, - ) - - self._component_collections = self.create_component_collections() + # create_component_collections creates the per-Q parameter lists (A_0/A_1 and + # lorentzian_width) itself, so the components it builds are backed by the very + # parameters stored in those lists. + self.create_component_collections() # ------------------------------------------------------------------ # Properties @@ -492,7 +483,13 @@ def create_component_collections( self, ) -> list[ComponentCollection]: r""" - Create ComponentCollections for the DeltaLorentz model at given Q values. + Create ComponentCollections for the DeltaLorentz model at given Q values. + + The per-Q parameter lists (A_0/A_1 and lorentzian_width, when Q-variation is enabled) are + recreated here so the built components are backed by the very parameters stored in the + lists, keeping ``calculate_width``/``calculate_EISF``/``calculate_QISF`` in sync with the + components. The created collections are installed on the model (they become the collections + returned by ``get_component_collections``), so the returned list is the live one. Returns ------- @@ -501,24 +498,31 @@ def create_component_collections( value. """ if self.Q is None: - return [] + self._A_0_list = [] + self._A_1_list = [] + self._lorentzian_width_list = [] + self._component_collections = [] + return self._component_collections Q = self.Q.values if self._allow_Q_variation['A_0'] is True: - A_0_list, A_1_list = self._create_A0_A1_parameter_lists(self.A_0) - self._A_0_list = A_0_list - self._A_1_list = A_1_list + self._A_0_list, self._A_1_list = self._create_A0_A1_parameter_lists(self.A_0) + else: + self._A_0_list = [] + self._A_1_list = [] if self._allow_Q_variation['lorentzian_width'] is True: - lorentzian_width_list = self._create_lorentzian_width_parameter_list( + self._lorentzian_width_list = self._create_lorentzian_width_parameter_list( self.lorentzian_width ) - self._lorentzian_width_list = lorentzian_width_list + else: + self._lorentzian_width_list = [] component_collection_list = [None] * len(Q) for i, Q_value in enumerate(Q): component_collection_list[i] = ComponentCollection( + name=f'{self.name}_Q{Q_value:.2f}', display_name=f'{self.display_name}_Q{Q_value:.2f}', x_unit=self.x_unit, y_unit=self.y_unit, @@ -588,7 +592,8 @@ def create_component_collections( component_collection_list[i].append_component(delta_component) - return component_collection_list + self._component_collections = component_collection_list + return self._component_collections def get_fit_targets(self) -> list[FitTarget]: """ @@ -877,8 +882,13 @@ def _create_A0_A1_parameter_lists( A_0_list = [] A_1_list = [] for _ in range(len(self.Q)): + # Like the per-Q width parameters (named ' width'), the per-Q + # amplitudes carry the model name so they do not collide with other models' + # parameters. The name is the same at every Q on purpose: parameters are tracked + # across Q by name (unique within a Q, shared across Q). a0 = Parameter( - name='A_0', + name=f'{self.name} A_0', + display_name='A_0', value=float(A_0.value), fixed=False, min=0.0, @@ -886,7 +896,7 @@ def _create_A0_A1_parameter_lists( ) a1 = Parameter.from_dependency( - name='A_1', + name=f'{self.name} A_1', dependency_expression='1 - A_0', dependency_map={'A_0': a0}, ) @@ -931,27 +941,12 @@ def _create_lorentzian_width_parameter_list( def _on_Q_change(self) -> None: """ - Handle changes to the Q values. Updates the A_0, A_1 and lorentzian_width parameters if - they are allowed to vary with Q. - """ - if self.Q is None: - self._A_0_list = [] - self._A_1_list = [] - self._lorentzian_width_list = [] - else: - if self._allow_Q_variation['A_0'] is True: - self._A_0_list, self._A_1_list = self._create_A0_A1_parameter_lists(self.A_0) - else: - self._A_0_list = [] - self._A_1_list = [] + Handle changes to the Q values. - if self._allow_Q_variation['lorentzian_width'] is True: - self._lorentzian_width_list = self._create_lorentzian_width_parameter_list( - self.lorentzian_width - ) - else: - self._lorentzian_width_list = [] - self._component_collections = self.create_component_collections() + Rebuilds the component collections; the per-Q A_0, A_1 and lorentzian_width parameter lists + are recreated inside ``create_component_collections``. + """ + self.create_component_collections() def _convert_extra_x_unit_parameters(self, unit_str: str) -> None: """ diff --git a/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py b/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py index c8bafd083..0e4b3ed2b 100644 --- a/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py +++ b/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py @@ -203,8 +203,11 @@ def Q(self, value: Q_type | None) -> None: if len(old_Q) != len(new_Q) or not sc.allclose(old_Q, new_Q): raise ValueError( - 'New Q values are not similar to the old ones. ' - 'To change Q values, first run clear_Q().' + f'New Q values are not similar to the old ones on diffusion model ' + f'{self.name!r}. This typically happens when a diffusion model that was ' + f'previously used with different Q values (e.g. in another SampleModel) is ' + f'reused. Run clear_Q(confirm=True) on the diffusion model first, then set ' + f'the new Q values.' ) @property diff --git a/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py b/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py index ad615e3b2..5c29ceb73 100644 --- a/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py +++ b/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py @@ -22,7 +22,7 @@ class JumpTranslationalDiffusion(DiffusionModelBase): The model consists of a Lorentzian function for each Q-value, where the width is given by - $$ \Gamma(Q) = \frac{Q^2}{1+D t Q^2}. $$ + $$ \Gamma(Q) = \frac{\hbar D Q^2}{1+D t Q^2}. $$ where $D$ is the diffusion coefficient and $t$ is the relaxation time. Q is assumed to have units of 1/angstrom. Creates ComponentCollections with Lorentzian components for given @@ -87,10 +87,10 @@ def __init__( Display name of the diffusion model. lorentzian_name : str | None, default=None Name of the Lorentzian component. If None, it will be set to the name of the diffusion - model with '_Lorentzian' appended. By default, None. + model. By default, None. lorentzian_display_name : str | None, default=None - Display name of the Lorentzian component. If None, it will be set to the display name - of the diffusion model with '_Lorentzian' appended. By default, None + Display name of the Lorentzian component. If None, it will be set to the + lorentzian_name. By default, None unique_name : str | None, default=None Unique name of the diffusion model. If None, a unique name will be generated. By default, None. @@ -305,6 +305,9 @@ def create_component_collections( """ Create ComponentCollection components for the diffusion model at given Q values. + The created collections are installed on the model (they become the collections returned by + ``get_component_collections``), so the returned list is the live one. + Returns ------- list[ComponentCollection] @@ -361,7 +364,8 @@ def create_component_collections( component_collection_list[i].append_component(lorentzian_component) - return component_collection_list + self._component_collections = component_collection_list + return self._component_collections ################################ # Private methods diff --git a/src/easydynamics/sample_model/model_base.py b/src/easydynamics/sample_model/model_base.py index 0e0e51cc7..15040a852 100644 --- a/src/easydynamics/sample_model/model_base.py +++ b/src/easydynamics/sample_model/model_base.py @@ -75,6 +75,13 @@ def __init__( self._components = ComponentCollection(x_unit=self.x_unit, y_unit=self.y_unit) self._component_collections: list[ComponentCollection] = [] + # Counter part of state_version: bumped whenever the dirty flag is raised. + self._state_counter = 0 + # Template-collection version the per-Q collections were last built from. Compared + # against self._components.version so in-place mutations of the live template + # collection (reachable via the `components` property) are detected without + # callbacks. + self._built_components_version = self._components.version self._component_collections_is_dirty = True if isinstance(components, (ModelComponent, ComponentCollection)): self.append_component(components) @@ -97,7 +104,7 @@ def evaluate( Raises ------ ValueError - If there are no components in the model to evaluate. + If Q is not set on the model, or if there are no components in the model to evaluate. Returns ------- @@ -107,6 +114,11 @@ def evaluate( """ self._ensure_component_collections_current() if not self._component_collections: + if self.Q is None: + raise ValueError( + 'Q is not set on the model, so there are no per-Q component collections ' + 'to evaluate. Set Q before evaluating.' + ) raise ValueError('No components in the model to evaluate.') return [ collection.evaluate(x, output=output) for collection in self._component_collections @@ -149,14 +161,18 @@ def clear_components(self) -> None: # ------------------------------------------------------------------ @property - def components(self) -> list[ModelComponent]: + def components(self) -> ComponentCollection: """ - Get the components of the SampleModel. + Get the template ComponentCollection of the SampleModel. + + This is the live template collection: mutating it in place (e.g. via ``append_component``) + is detected through its ``version`` and triggers a rebuild of the per-Q collections on next + use. Returns ------- - list[ModelComponent] - The components of the SampleModel. + ComponentCollection + The template component collection of the SampleModel. """ return self._components @@ -187,12 +203,72 @@ def component_collections_is_dirty(self) -> bool: """ Return whether component collections need to be rebuilt before use. + Collections are stale when the dirty flag was raised (Q or component changes through the + model's methods) or when the live template collection was mutated in place since the + collections were last built. + Returns ------- bool ``True`` if component collections have not been built yet or are stale. """ - return self._component_collections_is_dirty + return ( + self._component_collections_is_dirty + or self._built_components_version != self._components.version + ) + + @property + def _component_collections_is_dirty(self) -> bool: + """ + Get the dirty flag for the per-Q component collections. + + Implemented as a property so every write is intercepted: raising the flag bumps the state + counter (making ``state_version`` change), and clearing it records the template collection + version the collections were built from. + + Returns + ------- + bool + The raw dirty flag (does not account for in-place template mutations; use + ``component_collections_is_dirty`` for the full staleness check). + """ + return self._component_collections_dirty_flag + + @_component_collections_is_dirty.setter + def _component_collections_is_dirty(self, value: bool) -> None: + """ + Set the dirty flag for the per-Q component collections. + + Parameters + ---------- + value : bool + ``True`` marks the collections stale and bumps the state counter. ``False`` marks them + current and records the template collection version they now correspond to. + """ + value = bool(value) + if value: + self._state_counter += 1 + else: + self._built_components_version = self._components.version + self._component_collections_dirty_flag = value + + @property + def state_version(self) -> int: + """ + Get a monotonic version of everything affecting the per-Q component collections. + + The value changes whenever Q changes, components are added/removed/replaced through the + model's methods, or the live template collection (``components``) is mutated in place. + Implemented as an internal counter plus the template collection's mutation version, so it + only ever increases. Reading never rebuilds, clears or mutates anything; equal values mean + the collections' inputs are unchanged. + + Returns + ------- + int + The current state version. + """ + return self._state_counter + self._components.version @property def Q(self) -> sc.Variable | None: @@ -396,8 +472,11 @@ def normalize_area(self) -> None: def _ensure_component_collections_current(self) -> None: """ Rebuild component collections if any dependency has changed since they were last built. + + Uses the full staleness check, so both flag-raising changes (Q, component methods) and + in-place mutations of the live template collection trigger a rebuild. """ - if self._component_collections_is_dirty: + if self.component_collections_is_dirty: self._generate_component_collections() self._component_collections_is_dirty = False diff --git a/src/easydynamics/sample_model/resolution_model.py b/src/easydynamics/sample_model/resolution_model.py index 3949f843c..b33024ab3 100644 --- a/src/easydynamics/sample_model/resolution_model.py +++ b/src/easydynamics/sample_model/resolution_model.py @@ -1,6 +1,7 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import warnings from copy import copy import scipp as sc @@ -74,6 +75,8 @@ def __init__( Q : Q_type | None, default=None Q values for the model. If None, Q is not set. """ + # Set before super().__init__, which may call append_component (overridden below). + self._calibrated = False super().__init__( display_name=display_name, unique_name=unique_name, @@ -99,7 +102,14 @@ def append_component(self, component: ModelComponent | ComponentCollection) -> N ------ TypeError If the component is a DeltaFunction, Polynomial, or Exponential. + + Notes + ----- + A ``RuntimeError`` propagates from the calibration guard if the model holds calibrated + per-Q collections from ``from_sample_model``; a template change would schedule a rebuild + that silently discards them. """ + self._assert_not_calibrated('append a component') components = component if isinstance(component, ComponentCollection) else (component,) for comp in components: @@ -110,6 +120,84 @@ def append_component(self, component: ModelComponent | ComponentCollection) -> N super().append_component(component) + def remove_component(self, name: str) -> None: + """ + Remove a component from the ResolutionModel by its name. + + Parameters + ---------- + name : str + The name of the component to remove. + + Notes + ----- + A ``RuntimeError`` propagates from the calibration guard if the model holds calibrated + per-Q collections from ``from_sample_model``; a template change would schedule a rebuild + that silently discards them. + """ + self._assert_not_calibrated('remove a component') + super().remove_component(name) + + def clear_components(self) -> None: + """ + Clear all components from the ResolutionModel. + + Notes + ----- + A ``RuntimeError`` propagates from the calibration guard if the model holds calibrated + per-Q collections from ``from_sample_model``; a template change would schedule a rebuild + that silently discards them. + """ + self._assert_not_calibrated('clear the components') + super().clear_components() + + def clear_Q(self, confirm: bool = False) -> None: + """ + Clear the Q values of the ResolutionModel, removing all component collections and their + associated Parameters. + + Parameters + ---------- + confirm : bool, default=False + Confirmation to clear Q values. + + Notes + ----- + A ``ValueError`` propagates from the base implementation if confirm is not True, and a + ``RuntimeError`` propagates from the calibration guard if the model holds calibrated per-Q + collections from ``from_sample_model``; clearing Q would discard them. + """ + self._assert_not_calibrated('clear Q') + super().clear_Q(confirm=confirm) + + def _assert_not_calibrated(self, action: str) -> None: + """ + Raise if this model holds calibrated per-Q collections installed by from_sample_model. + + The per-Q collections installed by ``from_sample_model`` hold the fitted (calibrated) + resolution, but the template components do not. Any mutation that schedules a rebuild would + silently replace the calibrated collections with unfitted template copies, so such + mutations fail loudly instead. + + Parameters + ---------- + action : str + Description of the attempted mutation, used in the error message. + + Raises + ------ + RuntimeError + If the model is calibrated. + """ + if self._calibrated: + raise RuntimeError( + f'Cannot {action} on a ResolutionModel created by from_sample_model: its per-Q ' + f'collections hold the fitted (calibrated) resolution, and this change would ' + f'rebuild them from the unfitted template, silently discarding the calibration. ' + f'Create a new ResolutionModel (or rerun from_sample_model on an updated ' + f'SampleModel) instead.' + ) + @classmethod def from_sample_model( cls, @@ -120,6 +208,18 @@ def from_sample_model( """ Create a ResolutionModel from a SampleModel. + DeltaFunction components (the standard QENS elastic line) are stripped from both the + template and the per-Q collections, with a warning: a delta carries no resolution + broadening (it is the identity under convolution), so the fitted broadened components are + the resolution. Polynomial and Exponential components are rejected, as backgrounds do not + belong in a resolution model. + + When the SampleModel has Q values, the fitted per-Q collections are installed as the + calibrated resolution and the model is locked: mutations that would rebuild the collections + from the (unfitted) template — ``append_component``, ``remove_component``, + ``clear_components``, ``clear_Q`` — raise a RuntimeError instead of silently discarding the + calibration. + Parameters ---------- sample_model : SampleModel @@ -137,8 +237,11 @@ def from_sample_model( Raises ------ TypeError - If sample_model is not a SampleModel, or if normalize_area or fix_parameters are not - bool. + If sample_model is not a SampleModel, if normalize_area or fix_parameters are not bool, + or if the SampleModel contains Polynomial or Exponential components. + ValueError + If a per-Q collection contains only DeltaFunction components, leaving no resolution + shape after stripping. """ if not isinstance(sample_model, SampleModel): raise TypeError( @@ -151,11 +254,22 @@ def from_sample_model( if not isinstance(fix_parameters, bool): raise TypeError('fix_parameters must be True or False.') + template = ComponentCollection( + x_unit=sample_model.x_unit, + y_unit=sample_model.y_unit, + ) + stripped_deltas = 0 + for component in sample_model.components: + if isinstance(component, DeltaFunction): + stripped_deltas += 1 + continue + template.append_component(component) + resolution_model = cls( display_name=sample_model.display_name, x_unit=sample_model.x_unit, y_unit=sample_model.y_unit, - components=sample_model.components, + components=template, Q=sample_model.Q, ) @@ -163,10 +277,33 @@ def from_sample_model( # Prepare the per-Q collections detached from the model so no EasyScience # callback can schedule a rebuild halfway through, then install them and # clear the dirty flag in one final step. - collections = [ - copy(sample_model.get_component_collection(Q_index=index)) - for index in range(len(sample_model.Q)) - ] + collections = [] + for index in range(len(sample_model.Q)): + source = copy(sample_model.get_component_collection(Q_index=index)) + filtered = ComponentCollection( + name=source.name, + display_name=source.display_name, + x_unit=source.x_unit, + y_unit=source.y_unit, + ) + for component in source: + if isinstance(component, DeltaFunction): + stripped_deltas += 1 + continue + if isinstance(component, (Polynomial, Exponential)): + raise TypeError( + f'Component in ResolutionModel cannot be a ' + f'{component.__class__.__name__}' + ) + filtered.append_component(component) + if len(filtered) == 0: + raise ValueError( + f'The SampleModel collection at Q index {index} contains only ' + f'DeltaFunction components; after stripping them no resolution shape ' + f'is left. Fit the resolution data with at least one broadened ' + f'component (e.g. a Gaussian).' + ) + collections.append(filtered) for collection in collections: if normalize_area: collection.normalize_area() @@ -174,6 +311,16 @@ def from_sample_model( collection.fix_all_parameters() resolution_model._component_collections = collections resolution_model._component_collections_is_dirty = False + resolution_model._calibrated = True + + if stripped_deltas: + warnings.warn( + f'Stripped {stripped_deltas} DeltaFunction component(s) from the SampleModel ' + f'when building the ResolutionModel: a delta function carries no resolution ' + f'broadening (it is the identity under convolution).', + UserWarning, + stacklevel=2, + ) return resolution_model diff --git a/src/easydynamics/sample_model/sample_model.py b/src/easydynamics/sample_model/sample_model.py index fc6b77624..20ae75da7 100644 --- a/src/easydynamics/sample_model/sample_model.py +++ b/src/easydynamics/sample_model/sample_model.py @@ -16,6 +16,7 @@ from easydynamics.utils.utils import Numeric from easydynamics.utils.utils import Q_type from easydynamics.utils.utils import _validate_and_convert_Q +from easydynamics.utils.utils import _validate_unit from easydynamics.utils.utils import convert_units_with_rollback @@ -125,42 +126,48 @@ def __init__( self._diffusion_models = diffusion_models Q = _validate_and_convert_Q(Q) - for dm in self.diffusion_models: - dm.Q = Q # Ensure diffusion models have the same Q as the SampleModel - - super().__init__( - display_name=display_name, - unique_name=unique_name, - x_unit=x_unit, - y_unit=y_unit, - components=components, - Q=Q, - ) + # Validate (and build) everything else before mutating the passed diffusion models + # below, so a failed construction does not leave them changed (their Q set and their + # component collections rebuilt). + temperature_unit = _validate_unit(temperature_unit) if temperature is None: - self._temperature = None + temperature_parameter = None else: if not isinstance(temperature, Numeric): raise TypeError('temperature must be a number or None') if temperature < 0: raise ValueError('temperature must be non-negative') - self._temperature = Parameter( + temperature_parameter = Parameter( name='Temperature', value=temperature, unit=temperature_unit, display_name='Temperature', fixed=True, ) - self._temperature_unit = temperature_unit if detailed_balance_settings is None: - self._detailed_balance_settings = DetailedBalanceSettings() - elif isinstance(detailed_balance_settings, DetailedBalanceSettings): - self._detailed_balance_settings = detailed_balance_settings - else: + detailed_balance_settings = DetailedBalanceSettings() + elif not isinstance(detailed_balance_settings, DetailedBalanceSettings): raise TypeError('detailed_balance_settings must be a DetailedBalanceSettings or None') + for dm in self.diffusion_models: + dm.Q = Q # Ensure diffusion models have the same Q as the SampleModel + + super().__init__( + display_name=display_name, + unique_name=unique_name, + x_unit=x_unit, + y_unit=y_unit, + components=components, + Q=Q, + ) + + self._temperature = temperature_parameter + self._temperature_unit = temperature_unit + self._detailed_balance_settings = detailed_balance_settings + # ------------------------------------------------------------------ # Component management # ------------------------------------------------------------------ @@ -326,14 +333,14 @@ def temperature(self, value: Numeric | None) -> None: self._temperature.value = value @property - def temperature_unit(self) -> str | sc.Unit: + def temperature_unit(self) -> str: """ Get the temperature unit. Returns ------- - str | sc.Unit - The unit of the temperature parameter. + str + The unit of the temperature parameter, normalized to a string. """ return self._temperature_unit @@ -378,6 +385,7 @@ def convert_temperature_unit(self, unit: str | sc.Unit) -> None: if self.temperature is None: raise ValueError('Temperature is not set, cannot convert unit.') + unit = _validate_unit(unit) # normalize to str, as easyscience expects old_unit = self.temperature.unit try: @@ -555,7 +563,13 @@ def evaluate( divide_by_temperature=self.detailed_balance_settings.normalize_detailed_balance, energy_unit=self.x_unit, ) - y = [yi * DBF for yi in y] + if output == 'scipp': + # DBF is a plain numpy array (a dimensionless factor when + # normalize_detailed_balance is True), so multiply the values and keep the + # unit label the collections produced, consistent with numpy output. + y = [sc.array(dims=yi.dims, values=yi.values * DBF, unit=yi.unit) for yi in y] + else: + y = [yi * DBF for yi in y] return y diff --git a/src/easydynamics/settings/convolution_settings.py b/src/easydynamics/settings/convolution_settings.py index 1ac059c73..46bdcb6ec 100644 --- a/src/easydynamics/settings/convolution_settings.py +++ b/src/easydynamics/settings/convolution_settings.py @@ -140,44 +140,51 @@ def upsample_factor(self, factor: Numeric | None) -> None: self._invalidate_plan() @property - def extension_factor(self) -> float: + def extension_factor(self) -> float | None: """ Get the extension factor. The extension factor determines how much the energy range is extended on both sides before - convolution. 0.2 means extending by 20% of the original energy span on each side + convolution. 0.2 means extending by 20% of the original energy span on each side. None is + only valid while upsampling is disabled (upsample_factor=None). Returns ------- - float - The extension factor. + float | None + The extension factor, or None if unset. """ return self._extension_factor @extension_factor.setter - def extension_factor(self, factor: Numeric) -> None: + def extension_factor(self, factor: Numeric | None) -> None: """ Set the extension factor and recreate the dense grid. The extension factor determines how much the energy range is extended on both sides before - convolution. 0.2 means extending by 20% of the original energy span on each side. + convolution. 0.2 means extending by 20% of the original energy span on each side. None is + accepted (matching the constructor), but convolvers require a numeric extension factor + whenever upsample_factor is set. Parameters ---------- - factor : Numeric + factor : Numeric | None The new extension factor. Raises ------ TypeError - If factor is not a number. + If factor is neither a number nor None. ValueError If factor is negative. """ + if factor is None: + self._extension_factor = None + self._invalidate_plan() + return if not isinstance(factor, Numeric): - raise TypeError('Extension factor must be a number.') + raise TypeError('Extension factor must be a number or None.') if factor < 0.0: raise ValueError('Extension factor must be non-negative.') diff --git a/src/easydynamics/settings/detailed_balance_settings.py b/src/easydynamics/settings/detailed_balance_settings.py index d8c4f7e9e..ea85c6e35 100644 --- a/src/easydynamics/settings/detailed_balance_settings.py +++ b/src/easydynamics/settings/detailed_balance_settings.py @@ -76,6 +76,41 @@ def __init__( unique_name=unique_name, ) + # Plan-invalidation bookkeeping for convolvers sharing this settings object. + # Mirrors ConvolutionSettings: _plan_version is bumped whenever a flag changes; + # each convolver records the version it last rebuilt against and rebuilds when the + # versions differ. + self._plan_version = 0 + + # ------------------------------------------------------------------ + # Plan invalidation + # ------------------------------------------------------------------ + + def _invalidate_plan(self) -> None: + """ + Invalidate the convolution plan for every convolver sharing these settings. + + Bumps the plan version, so every convolver that recorded an earlier version rebuilds its + plan before the next convolution. + """ + self._plan_version += 1 + + def _plan_valid_for(self, seen_version: int) -> bool: + """ + Check whether a convolver that last rebuilt at seen_version can skip rebuilding. + + Parameters + ---------- + seen_version : int + The plan version the convolver recorded when it last rebuilt its plan. + + Returns + ------- + bool + True if no invalidation happened since the convolver's rebuild. + """ + return seen_version == self._plan_version + # ------------------------------------------------------------------ # Properties # ------------------------------------------------------------------ @@ -110,6 +145,7 @@ def use_detailed_balance(self, value: bool) -> None: if not isinstance(value, bool): raise TypeError('use_detailed_balance must be True or False') self._use_detailed_balance = value + self._invalidate_plan() @property def normalize_detailed_balance(self) -> bool: @@ -143,6 +179,7 @@ def normalize_detailed_balance(self, value: bool) -> None: if not isinstance(value, bool): raise TypeError('normalize_detailed_balance must be True or False') self._normalize_detailed_balance = value + self._invalidate_plan() def __repr__(self) -> str: """ diff --git a/src/easydynamics/utils/__init__.py b/src/easydynamics/utils/__init__.py index 1c3402ced..bb782a682 100644 --- a/src/easydynamics/utils/__init__.py +++ b/src/easydynamics/utils/__init__.py @@ -6,9 +6,11 @@ from easydynamics.utils.posterior_plotting import plot_corner from easydynamics.utils.posterior_plotting import plot_posterior_predictive from easydynamics.utils.posterior_plotting import plot_trace +from easydynamics.utils.utils import hbar __all__ = [ 'detailed_balance_factor', + 'hbar', 'plot_corner', 'plot_posterior_predictive', 'plot_trace', diff --git a/src/easydynamics/utils/detailed_balance.py b/src/easydynamics/utils/detailed_balance.py index a3b540cdc..374f839a6 100644 --- a/src/easydynamics/utils/detailed_balance.py +++ b/src/easydynamics/utils/detailed_balance.py @@ -23,7 +23,7 @@ def detailed_balance_factor( - energy: float | list | np.ndarray | sc.Variable, + energy: float | list | np.ndarray | sc.Variable | sc.DataArray, temperature: float | sc.Variable | Parameter, energy_unit: str | sc.Unit = 'meV', temperature_unit: str | sc.Unit = 'K', @@ -37,10 +37,12 @@ def detailed_balance_factor( Parameters ---------- - energy : float | list | np.ndarray | sc.Variable - The energy transfer. If number, assumed to be in meV unless energy_unit is set. + energy : float | list | np.ndarray | sc.Variable | sc.DataArray + The energy transfer. If number, assumed to be in meV unless energy_unit is set. If a + DataArray, its single coordinate is used as the energy axis. temperature : float | sc.Variable | Parameter - The temperature. If number, assumed to be in K unless temperature_unit is set. + The temperature. Must be a single scalar value. If number, assumed to be in K unless + temperature_unit is set. energy_unit : str | sc.Unit, default='meV' Unit for energy if energy is given as a number or list. temperature_unit : str | sc.Unit, default='K' @@ -52,13 +54,13 @@ def detailed_balance_factor( Raises ------ TypeError - If energy or temperature is not a number, list, numpy array, or scipp Variable, or if - energy_unit or temperature_unit is not a string or scipp Unit, or if divide_by_temperature - is not a boolean. + If energy or temperature is not one of the accepted types, or if energy_unit or + temperature_unit is not a string or scipp Unit, or if divide_by_temperature is not a + boolean. ValueError - If temperature is negative, or if energy is a numpy array with more than 1 dimension, or if - temperature is a scipp Variable that does not have a single dimension named 'temperature', - or if energy is a scipp Variable that does not have a single dimension named 'energy'. + If temperature is negative or is not a single scalar value, if energy is a list or numpy + array with more than 1 dimension, or if energy is a scipp DataArray without exactly one + coordinate. UnitError If the provided energy_unit or temperature_unit is invalid, or if the units of energy or temperature cannot be converted to the expected units. @@ -109,6 +111,12 @@ def detailed_balance_factor( value=temperature, unit=temperature_unit, name='temperature' ) + if temperature.sizes != {}: + raise ValueError( + f'temperature must be a single scalar value, ' + f'got an array with sizes {dict(temperature.sizes)}.' + ) + if temperature.value < 0: raise ValueError('Temperature must be non-negative.') @@ -190,7 +198,7 @@ def detailed_balance_factor( def _convert_to_scipp_variable( - value: float | list | np.ndarray | Parameter | sc.Variable, + value: float | list | np.ndarray | Parameter | sc.Variable | sc.DataArray, name: str, unit: str | None = None, ) -> sc.Variable: @@ -199,9 +207,11 @@ def _convert_to_scipp_variable( Parameters ---------- - value : float | list | np.ndarray | Parameter | sc.Variable - The value to convert. Can be a number, list, numpy array, Parameter, or scipp Variable. If - a number or list, the unit must be specified in the unit argument. + value : float | list | np.ndarray | Parameter | sc.Variable | sc.DataArray + The value to convert. Can be a number, list, numpy array, Parameter, scipp Variable, or + scipp DataArray. If a number or list, the unit must be specified in the unit argument. A + DataArray must have exactly one coordinate, which is used as the value (consistent with how + components treat DataArray input to ``evaluate``). name : str The name of the variable, used for error messages. unit : str | None, default=None @@ -213,6 +223,9 @@ def _convert_to_scipp_variable( ------ TypeError If value is not one of the accepted types, or if unit is not a string when needed. + ValueError + If value is a list or numpy array with more than 1 dimension, or a DataArray without + exactly one coordinate. UnitError If the provided unit is invalid. @@ -221,6 +234,16 @@ def _convert_to_scipp_variable( sc.Variable The input value converted to a scipp Variable with appropriate units. """ + if isinstance(value, sc.DataArray): + coords = dict(value.coords) + if len(coords) != 1: + coord_names = ', '.join(coords.keys()) + raise ValueError( + f'scipp.DataArray must have exactly one coordinate to be used as {name}. ' + f'Found {len(coords)} coordinates: {coord_names}.' + ) + value = next(iter(coords.values())) + if isinstance(value, sc.Variable): return value @@ -237,6 +260,11 @@ def _convert_to_scipp_variable( raise TypeError(f'{name} must be a number, list, numpy array or scipp Variable') raise TypeError(f'{name} must be a number, list, numpy array, Parameter or scipp Variable') + if array_value.ndim > 1: + raise ValueError( + f'{name} must be at most one-dimensional, got {array_value.ndim} dimensions.' + ) + # Create appropriate scipp variable based on shape if array_value.shape == () or (array_value.shape == (1,)): # Scalar or single-element array diff --git a/src/easydynamics/utils/utils.py b/src/easydynamics/utils/utils.py index 06661df25..4ebc4a304 100644 --- a/src/easydynamics/utils/utils.py +++ b/src/easydynamics/utils/utils.py @@ -44,14 +44,15 @@ def verify_Q_index(Q_index: int, Q: sc.Variable | None, allow_none: bool = False Raises ------ TypeError - If Q_index is not an int (or not an int or None when allow_none=True). + If Q_index is not an int (or not an int or None when allow_none=True). Booleans are + rejected explicitly, since ``True`` would otherwise silently mean index 1. IndexError If Q_index is negative, or out of range when Q is available. """ if allow_none and Q_index is None: return - if Q_index is None or not isinstance(Q_index, int): + if Q_index is None or isinstance(Q_index, bool) or not isinstance(Q_index, int): if allow_none: raise TypeError(f'Q_index must be an int or None, got {type(Q_index).__name__}') raise TypeError(f'Q_index must be an int, got {type(Q_index).__name__}') @@ -275,6 +276,8 @@ def _in_notebook() -> bool: True if in a Jupyter notebook, False otherwise. """ try: + # Imported here deliberately: IPython may be absent at runtime, and the except + # clause below turns that into the answer "not a notebook". from IPython import get_ipython shell = get_ipython().__class__.__name__ diff --git a/tests/functional/test_dummy.py b/tests/functional/test_dummy.py index b45b191fe..6a2a5cd51 100644 --- a/tests/functional/test_dummy.py +++ b/tests/functional/test_dummy.py @@ -1,8 +1,22 @@ # SPDX-FileCopyrightText: 2025-2026 EasyDynamics contributors # SPDX-License-Identifier: BSD-3-Clause +import numpy as np -def test_dummy(): - calculated = 2 + 2 - expected = 4 - assert calculated == expected +import easydynamics as edyn + + +def test_smoke_build_and_evaluate_model(): + # WHEN a minimal sample model with a single Lorentzian component + lorentzian = edyn.Lorentzian(area=1.0, width=0.1) + model = edyn.SampleModel(components=lorentzian) + + # THEN evaluating the component on a small energy grid + energy = np.linspace(-1.0, 1.0, 101) + y = lorentzian.evaluate(energy) + + # EXPECT the package installs, the model builds, and the evaluation is finite and peaked + assert model is not None + assert y.shape == energy.shape + assert np.all(np.isfinite(y)) + assert y.max() > 0.0 diff --git a/tests/integration/fitting/test_bayesian_sampling.py b/tests/integration/fitting/test_bayesian_sampling.py index 58fc48ab6..81c98ea57 100644 --- a/tests/integration/fitting/test_bayesian_sampling.py +++ b/tests/integration/fitting/test_bayesian_sampling.py @@ -18,6 +18,8 @@ mpl.use('Agg') +import matplotlib.pyplot as plt + from easydynamics.analysis.analysis1d import Analysis1d from easydynamics.experiment import Experiment from easydynamics.sample_model import InstrumentModel @@ -192,8 +194,6 @@ def test_subset_sampling_produces_a_single_column(self): def test_plots_render(self, sampled_analysis): # WHEN - import matplotlib.pyplot as plt - n_parameters = len(sampled_analysis.get_free_parameters()) # THEN diff --git a/tests/integration/fitting/test_bayesian_sampling_multi_q.py b/tests/integration/fitting/test_bayesian_sampling_multi_q.py index 74061b57e..52a4a93d8 100644 --- a/tests/integration/fitting/test_bayesian_sampling_multi_q.py +++ b/tests/integration/fitting/test_bayesian_sampling_multi_q.py @@ -19,6 +19,8 @@ mpl.use('Agg') +import matplotlib.pyplot as plt + import easydynamics as edyn import easydynamics.sample_model as sm @@ -142,8 +144,6 @@ def test_sampling_leaves_the_fitted_values_untouched(self): def test_plots_render(self, simultaneously_sampled): # WHEN - import matplotlib.pyplot as plt - n_parameters = len(simultaneously_sampled._chain_parameters()) # THEN diff --git a/tests/integration/fitting/test_fitting_with_diffusion_model.py b/tests/integration/fitting/test_fitting_with_diffusion_model.py index f9912c249..067692de1 100644 --- a/tests/integration/fitting/test_fitting_with_diffusion_model.py +++ b/tests/integration/fitting/test_fitting_with_diffusion_model.py @@ -3,6 +3,7 @@ import numpy as np import pooch +import pytest from easydynamics.analysis.analysis import Analysis from easydynamics.experiment import Experiment @@ -19,6 +20,10 @@ from easydynamics.sample_model.resolution_model import ResolutionModel from easydynamics.sample_model.sample_model import SampleModel +# Every test here downloads its data files through pooch; deselect with -m 'not network' +# when running offline. +pytestmark = pytest.mark.network + class TestFittingWithDiffusionModel: def test_fitting_with_diffusion_model(self): @@ -146,7 +151,7 @@ def test_fitting_with_diffusion_model(self): pars = diffusion_model.get_all_parameters() - tol = 10 * pars[0].error + tol = 3 * pars[0].error assert np.isclose(pars[0].value, 1.1258025622851164e-08, atol=tol) - tol = 10 * pars[1].error + tol = 3 * pars[1].error assert np.isclose(pars[1].value, 0.6937774083152299, atol=tol) diff --git a/tests/unit/easydynamics/analysis/test_analysis.py b/tests/unit/easydynamics/analysis/test_analysis.py index b8bda4f69..ab34bf879 100644 --- a/tests/unit/easydynamics/analysis/test_analysis.py +++ b/tests/unit/easydynamics/analysis/test_analysis.py @@ -1,12 +1,15 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +from collections import Counter from unittest.mock import MagicMock from unittest.mock import patch import numpy as np import pytest import scipp as sc +from easyscience.fitting.multi_fitter import MultiFitter +from easyscience.variable import Parameter import easydynamics as edyn import easydynamics.sample_model as sm @@ -16,6 +19,7 @@ from easydynamics.sample_model import SampleModel from easydynamics.sample_model.components.gaussian import Gaussian from easydynamics.settings.convolution_settings import ConvolutionSettings +from easydynamics.settings.detailed_balance_settings import DetailedBalanceSettings Q_VALUES = [0.5, 1.0, 1.5] @@ -466,11 +470,19 @@ def test_plot_data_and_model_with_residuals( @pytest.mark.parametrize('include_residuals', [True, False]) def test_data_and_model_to_datagroup(self, analysis, include_residuals): - # WHEN + # WHEN a custom energy grid is passed energy = sc.array(dims=['energy'], values=[20.0, 30.0, 40.0], unit='meV') - datagroup = analysis.data_and_model_to_datagroup( - energy=energy, include_residuals=include_residuals - ) + + # THEN residuals cannot be computed on a custom grid, so they are omitted with a warning + if include_residuals: + with pytest.warns(UserWarning, match='omitted'): + datagroup = analysis.data_and_model_to_datagroup( + energy=energy, include_residuals=include_residuals + ) + else: + datagroup = analysis.data_and_model_to_datagroup( + energy=energy, include_residuals=include_residuals + ) # EXPECT assert isinstance(datagroup, sc.DataGroup) @@ -478,12 +490,20 @@ def test_data_and_model_to_datagroup(self, analysis, include_residuals): assert 'Model' in datagroup assert sc.identical(datagroup['Data'], analysis.experiment.binned_data) assert sc.identical(datagroup['Model'], analysis._create_model_array(energy=energy)) - if include_residuals: - assert 'Residuals' in datagroup - assert sc.identical( - datagroup['Residuals'], - analysis.experiment.binned_data - analysis._create_model_array(), - ) + assert 'Residuals' not in datagroup + + def test_data_and_model_to_datagroup_residuals_on_experiment_grid(self, analysis): + # WHEN no custom energy grid is given + + # THEN + datagroup = analysis.data_and_model_to_datagroup(include_residuals=True) + + # EXPECT residuals present and consistent with the data and model on the same grid + assert 'Residuals' in datagroup + assert sc.identical( + datagroup['Residuals'], + analysis.experiment.binned_data - analysis._create_model_array(), + ) def test_data_and_model_to_datagroup_no_data_raises(self, analysis): # WHEN @@ -800,6 +820,137 @@ def test_on_convolution_settings_changed(self, analysis): assert analysis1d.convolution_settings.upsample_factor == 7 assert analysis1d.convolution_settings.extension_factor == pytest.approx(0.3) + def test_on_detailed_balance_settings_changed(self, analysis): + # WHEN the analysis list has been built with the old settings + _ = analysis.analysis_list + assert analysis._analysis_list_is_dirty is False + new_settings = DetailedBalanceSettings( + use_detailed_balance=False, normalize_detailed_balance=False + ) + + # THEN (this calls _on_detailed_balance_settings_changed internally) + analysis.detailed_balance_settings = new_settings + + # EXPECT the parent holds the new settings object and the per-Q analyses are rebuilt + # around it, so the change actually reaches every Q index + assert analysis.detailed_balance_settings is new_settings + assert analysis._analysis_list_is_dirty is True + for analysis1d in analysis.analysis_list: + assert analysis1d.detailed_balance_settings is new_settings + + def test_detailed_balance_settings_change_invalidates_the_fitter(self, analysis): + # WHEN + original = analysis.fitter + + # THEN + analysis.detailed_balance_settings = DetailedBalanceSettings(use_detailed_balance=False) + + # EXPECT + assert analysis.fitter is not original + + def test_rebin_invalidates_the_fitter_and_the_sampler(self, analysis): + # WHEN the fitter and sampler exist from before the rebin + original_fitter = analysis.fitter + sampler = analysis.bayesian + + # THEN - energy rebin leaves Q unchanged, so no confirm required + with ( + patch.object(analysis.experiment, 'rebin'), + patch.object(sampler, 'invalidate') as mock_invalidate, + ): + analysis.rebin({'energy': 2}) + + # EXPECT neither keeps referencing the pre-rebin Analysis1d objects and data + assert analysis.fitter is not original_fitter + mock_invalidate.assert_called_once() + + def test_simultaneous_fit_uses_the_configured_fitter(self, analysis): + # WHEN the cached fitter has been configured (e.g. its minimizer switched) + fake_fitter = MagicMock() + fake_fitter.fit.return_value = 'simultaneous_result' + analysis._fitter = fake_fitter + analysis._fitter_is_dirty = False + + # THEN + result = analysis.fit(fit_method='simultaneous') + + # EXPECT the configured fitter object performed the fit, not a throwaway MultiFitter + fake_fitter.fit.assert_called_once() + assert result == 'simultaneous_result' + + def test_uses_a_multifitter(self, multi_q_analysis): + # EXPECT + assert isinstance(multi_q_analysis.fitter, MultiFitter) + assert len(multi_q_analysis.fitter.fit_object) == len(Q_VALUES) + + def test_get_all_variables(self, analysis): + # WHEN + extra_par = Parameter(name='extra_par', value=1.0) + analysis._extra_parameters = [extra_par] + + # THEN + variables = analysis.get_all_variables() + + # EXPECT variables across every Q index plus the extra parameters + expected = analysis.sample_model.get_all_variables() + expected.extend(analysis.instrument_model.get_all_variables()) + expected.append(extra_par) + assert Counter(variables) == Counter(expected) + + def test_get_all_variables_on_an_empty_analysis(self): + # WHEN + analysis = Analysis(display_name='Empty') + + # THEN EXPECT no failure and no variables + assert analysis.get_all_variables() == [] + assert analysis.get_parameters_near_bounds() == [] + + def test_get_parameters_near_bounds_builds_no_fitter_or_sampler(self, analysis): + # WHEN neither the fitter nor the sampler exists yet + assert analysis._fitter is None + assert analysis._bayesian is None + + # THEN + analysis.get_parameters_near_bounds() + + # EXPECT listing parameters did not build them as side effects + assert analysis._fitter is None + assert analysis._bayesian is None + + ############# + # The bayesian sampler (the Analysis side of the contract) + ############# + + def test_bayesian_returns_the_cached_sampler(self, analysis): + # THEN + sampler = analysis.bayesian + + # EXPECT the same object on second access + assert sampler is analysis.bayesian + + def test_bayesian_is_invalidated_when_the_experiment_changes(self, analysis): + # WHEN + sampler = analysis.bayesian + new_experiment = Experiment(data=analysis.experiment.data.copy(deep=True)) + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis.experiment = new_experiment + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_when_the_sample_model_changes(self, analysis): + # WHEN + sampler = analysis.bayesian + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis.sample_model = SampleModel(components=Gaussian()) + + # EXPECT + mock_invalidate.assert_called_once() + def test_fit_single_Q_valid(self, analysis): # WHEN analysis.analysis_list[1].fit = MagicMock(return_value='fit_result_Q1') @@ -1173,6 +1324,13 @@ def test_repr(self, analysis): assert 'display_name=' in repr_str assert 'n_analyses=' in repr_str + def test_repr_reports_a_current_analysis_count(self, analysis): + # WHEN the analysis list has not been built yet + assert analysis._analysis_list == [] + + # THEN EXPECT repr ensures the list is current instead of reporting a stale count + assert 'n_analyses=3' in repr(analysis) + ############# # Chain parameters and labels ############# @@ -1243,8 +1401,6 @@ def test_single_q_analysis_keeps_plain_names(self): def test_parameter_from_outside_the_analysis_keeps_its_name(self, multi_q_analysis): # WHEN a parameter belongs to no Q index of this analysis - from easyscience.variable import Parameter - stranger = Parameter(name='Gaussian width', value=1.0) # EXPECT it is returned unqualified rather than mislabelled diff --git a/tests/unit/easydynamics/analysis/test_analysis1d.py b/tests/unit/easydynamics/analysis/test_analysis1d.py index 99ab53ab5..40ae70a0b 100644 --- a/tests/unit/easydynamics/analysis/test_analysis1d.py +++ b/tests/unit/easydynamics/analysis/test_analysis1d.py @@ -19,6 +19,12 @@ from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.sample_model.components.gaussian import Gaussian from easydynamics.sample_model.components.polynomial import Polynomial +from easydynamics.settings.detailed_balance_settings import DetailedBalanceSettings + +# The per-consumer convolver staleness tracking relies on the ModelBase.state_version contract; +# until it lands, the conservative fallback rebuilds on every prepare, so 'no rebuild' tests +# cannot pass. +HAS_STATE_VERSION = hasattr(SampleModel, 'state_version') class TestAnalysis1d: @@ -413,10 +419,16 @@ def test_data_and_model_to_datagroup(self, analysis1d, include_residuals): # WHEN energy = sc.array(dims=['energy'], values=[20.0, 30.0, 40.0], unit='meV') - # THEN - datagroup = analysis1d.data_and_model_to_datagroup( - energy=energy, include_residuals=include_residuals - ) + # THEN residuals cannot be computed on a custom grid, so they are omitted with a warning + if include_residuals: + with pytest.warns(UserWarning, match='omitted'): + datagroup = analysis1d.data_and_model_to_datagroup( + energy=energy, include_residuals=include_residuals + ) + else: + datagroup = analysis1d.data_and_model_to_datagroup( + energy=energy, include_residuals=include_residuals + ) # EXPECT assert isinstance(datagroup, sc.DataGroup) @@ -427,14 +439,20 @@ def test_data_and_model_to_datagroup(self, analysis1d, include_residuals): analysis1d.experiment.binned_data['Q', analysis1d.Q_index], ) assert sc.identical(datagroup['Model'], analysis1d._create_model_array(energy=energy)) - if include_residuals: - assert 'Residuals' in datagroup - assert sc.identical( - datagroup['Residuals'], - datagroup['Data'] - analysis1d._create_model_array(), - ) - else: - assert 'Residuals' not in datagroup + assert 'Residuals' not in datagroup + + def test_data_and_model_to_datagroup_residuals_on_experiment_grid(self, analysis1d): + # WHEN no custom energy grid is given + + # THEN + datagroup = analysis1d.data_and_model_to_datagroup(include_residuals=True) + + # EXPECT residuals present and consistent with the data and model on the same grid + assert 'Residuals' in datagroup + assert sc.identical( + datagroup['Residuals'], + datagroup['Data'] - analysis1d._create_model_array(), + ) def test_data_and_model_to_datagroup_no_data_raises(self, analysis1d): # WHEN @@ -1054,15 +1072,13 @@ def test_fit_marks_convolver_dirty_when_sample_model_components_change(self, ana # EXPECT - convolver was rebuilt (_ensure_convolver_current called _create_convolver) analysis1d._create_convolver.assert_called_once() + @pytest.mark.skipif(not HAS_STATE_VERSION, reason='pending ModelBase.state_version contract') def test_fit_does_not_rebuild_convolver_when_nothing_changed(self, analysis1d): """fit() should not call _create_convolver if nothing has changed since last fit.""" - # WHEN - build convolver and clear all dirty flags + # WHEN - a first fit has built the convolver against the current model state analysis1d._create_convolver = MagicMock(return_value=None) - analysis1d._convolver_is_dirty = False - analysis1d.sample_model._component_collections_is_dirty = False - analysis1d.instrument_model.resolution_model._component_collections_is_dirty = False - # THEN - call fit() with nothing changed + # THEN - fit once to sync, then fit again with nothing changed with patch( 'easydynamics.analysis.analysis1d.EasyScienceFitter', return_value=MagicMock(fit=MagicMock(return_value=MagicMock())), @@ -1076,8 +1092,10 @@ def test_fit_does_not_rebuild_convolver_when_nothing_changed(self, analysis1d): ) ) analysis1d.fit() + analysis1d._create_convolver.reset_mock() + analysis1d.fit() - # EXPECT - _create_convolver was NOT called (convolver reused) + # EXPECT - _create_convolver was NOT called again (convolver reused) analysis1d._create_convolver.assert_not_called() def test_rebin_rebins_experiment(self, analysis1d): @@ -1164,7 +1182,161 @@ def test_fit_marks_convolver_dirty_when_resolution_model_components_change(self, # EXPECT analysis1d._create_convolver.assert_called_once() - # ───── Regression tests ───── + ############# + # Convolver staleness across analyses sharing a model + ############# + + @pytest.fixture + def sibling_analyses(self): + """Two Analysis1d objects sharing one SampleModel and one InstrumentModel.""" + Q = sc.array(dims=['Q'], values=[1.0, 2.0], unit='1/Angstrom') + energy = sc.linspace('energy', -5.0, 5.0, num=11, unit='meV') + values = np.ones((2, 11)) + data_array = sc.DataArray( + data=sc.array(dims=['Q', 'energy'], values=values, variances=values), + coords={'Q': Q, 'energy': energy}, + ) + experiment = Experiment(data=data_array) + sample_model = SampleModel(components=Gaussian()) + instrument_model = InstrumentModel( + resolution_model=ResolutionModel(components=Gaussian(width=0.5)) + ) + return [ + Analysis1d( + display_name=f'Sibling{q_index}', + experiment=experiment, + sample_model=sample_model, + instrument_model=instrument_model, + Q_index=q_index, + ) + for q_index in (0, 1) + ] + + def test_in_place_model_edit_rebuilds_the_convolvers_of_all_siblings(self, sibling_analyses): + """Regression: the first sibling to prepare must not consume the staleness signal.""" + # WHEN both siblings have built their convolvers against the shared model + first, second = sibling_analyses + first._prepare_for_sampling() + second._prepare_for_sampling() + first_convolver = first._convolver + second_convolver = second._convolver + assert first_convolver is not None + assert second_convolver is not None + + # THEN the shared model is edited in place (not through any Analysis1d setter) + first.sample_model.append_component(Gaussian(name='ExtraGaussian')) + first._prepare_for_sampling() + second._prepare_for_sampling() + + # EXPECT both siblings rebuilt their convolvers, not only the first one to prepare + assert first._convolver is not first_convolver + assert second._convolver is not second_convolver + + def test_in_place_resolution_edit_rebuilds_the_convolvers_of_all_siblings( + self, sibling_analyses + ): + # WHEN both siblings have built their convolvers against the shared resolution model + first, second = sibling_analyses + first._prepare_for_sampling() + second._prepare_for_sampling() + first_convolver = first._convolver + second_convolver = second._convolver + + # THEN the shared resolution model is edited in place + first.instrument_model.resolution_model.append_component(Gaussian(name='ExtraResolution')) + first._prepare_for_sampling() + second._prepare_for_sampling() + + # EXPECT both siblings rebuilt their convolvers + assert first._convolver is not first_convolver + assert second._convolver is not second_convolver + + @pytest.mark.skipif(not HAS_STATE_VERSION, reason='pending ModelBase.state_version contract') + def test_prepare_does_not_rebuild_when_the_models_are_unchanged(self, sibling_analyses): + # WHEN a convolver has been built against the current model state + first, _ = sibling_analyses + first._prepare_for_sampling() + convolver = first._convolver + + # THEN preparing again with nothing changed + first._prepare_for_sampling() + + # EXPECT the convolver is reused, not rebuilt + assert first._convolver is convolver + + ############# + # Detailed balance settings + ############# + + def test_detailed_balance_settings_change_marks_convolver_dirty(self, analysis1d): + # WHEN + analysis1d._convolver_is_dirty = False + + # THEN a new settings object is assigned + analysis1d.detailed_balance_settings = DetailedBalanceSettings(use_detailed_balance=False) + + # EXPECT + assert analysis1d._convolver_is_dirty is True + + ############# + # The bayesian sampler (the Analysis1d side of the contract) + ############# + + def test_bayesian_returns_the_cached_sampler(self, analysis1d): + # THEN + sampler = analysis1d.bayesian + + # EXPECT the same object on second access + assert sampler is analysis1d.bayesian + + def test_bayesian_is_invalidated_when_the_Q_index_changes(self, analysis1d): + # WHEN + sampler = analysis1d.bayesian + + # THEN a different Q index means different data + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis1d.Q_index = 1 + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_when_the_experiment_changes(self, analysis1d): + # WHEN + sampler = analysis1d.bayesian + new_experiment = Experiment(data=analysis1d.experiment.data.copy(deep=True)) + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis1d.experiment = new_experiment + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_when_the_sample_model_changes(self, analysis1d): + # WHEN + sampler = analysis1d.bayesian + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis1d.sample_model = SampleModel(components=Gaussian()) + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_on_rebin(self, analysis1d): + # WHEN + sampler = analysis1d.bayesian + + # THEN rebinning changes the data the sampler was bound to + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis1d.rebin({'Q': 1}) + + # EXPECT + mock_invalidate.assert_called_once() + + ############# + # Regression tests + ############# @pytest.fixture def analysis1d_with_nan(self): @@ -1207,12 +1379,9 @@ def test_data_and_model_to_datagroup_with_nan_excludes_nan_from_data( # Before the fix, 'Data' contained the full 3-point grid (including NaN) # and computing Residuals crashed on the dimension mismatch. # WHEN - energy = sc.array(dims=['energy'], values=[20.0, 30.0, 40.0], unit='meV') # THEN - datagroup = analysis1d_with_nan.data_and_model_to_datagroup( - energy=energy, include_residuals=True - ) + datagroup = analysis1d_with_nan.data_and_model_to_datagroup(include_residuals=True) # EXPECT assert isinstance(datagroup, sc.DataGroup) @@ -1228,41 +1397,47 @@ def test_repr(self, analysis1d): assert 'display_name=' in repr_str assert 'Q_index=' in repr_str + ############# + # Change handlers + ############# -def _coverage_analysis1d(): - Q = sc.array(dims=['Q'], values=[1, 2, 3], unit='1/Angstrom') - energy = sc.array(dims=['energy'], values=[10.0, 20.0, 30.0], unit='meV') - data = sc.array( - dims=['Q', 'energy'], - values=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], - variances=[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]], - ) - data_array = sc.DataArray(data=data, coords={'Q': Q, 'energy': energy}) - return Analysis1d( - display_name='CoverageAnalysis', - experiment=Experiment(data=data_array), - sample_model=SampleModel(components=Gaussian()), - instrument_model=InstrumentModel(), - Q_index=0, - ) + @staticmethod + def _coverage_analysis1d(): + Q = sc.array(dims=['Q'], values=[1, 2, 3], unit='1/Angstrom') + energy = sc.array(dims=['energy'], values=[10.0, 20.0, 30.0], unit='meV') + data = sc.array( + dims=['Q', 'energy'], + values=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], + variances=[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]], + ) + data_array = sc.DataArray(data=data, coords={'Q': Q, 'energy': energy}) + return Analysis1d( + display_name='CoverageAnalysis', + experiment=Experiment(data=data_array), + sample_model=SampleModel(components=Gaussian()), + instrument_model=InstrumentModel(), + Q_index=0, + ) + + def test_on_Q_index_changed_with_none_clears_masked_energy(self): + # WHEN an analysis whose Q_index has been cleared + analysis1d = self._coverage_analysis1d() + analysis1d._Q_index = None + + # THEN the Q-index-changed handler runs + analysis1d._on_Q_index_changed() + + # EXPECT masked energy cleared and convolver marked dirty + assert analysis1d._masked_energy is None + assert analysis1d._convolver_is_dirty is True + + def test_on_experiment_changed_refreshes_masked_energy_when_Q_index_set(self): + # WHEN an analysis with a Q_index already set + analysis1d = self._coverage_analysis1d() + # THEN the experiment-changed handler runs + analysis1d._on_experiment_changed() -def test_on_Q_index_changed_with_none_clears_masked_energy(): - # GIVEN an analysis whose Q_index has been cleared - analysis1d = _coverage_analysis1d() - analysis1d._Q_index = None - # WHEN the Q-index-changed handler runs - analysis1d._on_Q_index_changed() - # EXPECT masked energy cleared and convolver marked dirty - assert analysis1d._masked_energy is None - assert analysis1d._convolver_is_dirty is True - - -def test_on_experiment_changed_refreshes_masked_energy_when_Q_index_set(): - # GIVEN an analysis with a Q_index already set - analysis1d = _coverage_analysis1d() - # WHEN the experiment-changed handler runs - analysis1d._on_experiment_changed() - # EXPECT masked energy refreshed and convolver marked dirty - assert analysis1d._masked_energy is not None - assert analysis1d._convolver_is_dirty is True + # EXPECT masked energy refreshed and convolver marked dirty + assert analysis1d._masked_energy is not None + assert analysis1d._convolver_is_dirty is True diff --git a/tests/unit/easydynamics/analysis/test_analysis_base.py b/tests/unit/easydynamics/analysis/test_analysis_base.py index afddabc83..499348934 100644 --- a/tests/unit/easydynamics/analysis/test_analysis_base.py +++ b/tests/unit/easydynamics/analysis/test_analysis_base.py @@ -378,6 +378,33 @@ def test_detailed_balance_settings_setter_invalid(self, analysis_base): ): analysis_base.detailed_balance_settings = 'invalid_settings' + def test_detailed_balance_settings_calls_on_detailed_balance_settings_changed( + self, analysis_base + ): + # WHEN + new_settings = DetailedBalanceSettings( + use_detailed_balance=False, normalize_detailed_balance=False + ) + with patch.object( + analysis_base, '_on_detailed_balance_settings_changed' + ) as mock_on_detailed_balance_settings_changed: + # THEN + analysis_base.detailed_balance_settings = new_settings + + # EXPECT the change hook fires, like every sibling setter's does + mock_on_detailed_balance_settings_changed.assert_called_once() + + def test_detailed_balance_settings_setter_invalid_fires_no_hook(self, analysis_base): + # WHEN / THEN + with ( + patch.object(analysis_base, '_on_detailed_balance_settings_changed') as mock_hook, + pytest.raises(TypeError), + ): + analysis_base.detailed_balance_settings = 'invalid_settings' + + # EXPECT + mock_hook.assert_not_called() + @pytest.mark.parametrize( 'extra_parameters', [ @@ -393,9 +420,6 @@ def test_detailed_balance_settings_setter_invalid(self, analysis_base): ], ) def test_extra_parameters_property(self, analysis_base, extra_parameters): - # WHEN - analysis_base.extra_parameters = extra_parameters - # THEN analysis_base.extra_parameters = extra_parameters diff --git a/tests/unit/easydynamics/analysis/test_fit_binding.py b/tests/unit/easydynamics/analysis/test_fit_binding.py index 901ab9e50..e1a4afcbe 100644 --- a/tests/unit/easydynamics/analysis/test_fit_binding.py +++ b/tests/unit/easydynamics/analysis/test_fit_binding.py @@ -26,9 +26,9 @@ def diffusion_binding(self): model = BrownianTranslationalDiffusion(lorentzian_name='Lorentzian') return FitBinding(model=model) - # ------------------------------------------------------------------ + ############# # Initialization and validation - # ------------------------------------------------------------------ + ############# def test_initialization(self, component_binding): # WHEN THEN EXPECT @@ -76,9 +76,9 @@ def test_diffusion_non_string_dataset_key_raises(self): with pytest.raises(TypeError, match='dataset keys'): FitBinding(model=model, targets={'width': 123}) - # ------------------------------------------------------------------ + ############# # Properties - # ------------------------------------------------------------------ + ############# def test_model_setter_revalidates_targets(self): # WHEN: a binding using DeltaLorentz's delta_area prediction @@ -105,9 +105,9 @@ def test_targets_setter_invalid_raises(self, diffusion_binding): with pytest.raises(ValueError, match='Unknown prediction'): diffusion_binding.targets = ['nonsense'] - # ------------------------------------------------------------------ + ############# # get_targets - # ------------------------------------------------------------------ + ############# def test_component_target(self, component_binding): # WHEN @@ -241,9 +241,9 @@ def test_delta_lorentz_delta_area_function(self): # EXPECT np.testing.assert_allclose(target.function(Q), model.calculate_EISF(Q) * model.scale.value) - # ------------------------------------------------------------------ + ############# # dunder methods - # ------------------------------------------------------------------ + ############# def test_repr(self, diffusion_binding): # WHEN THEN @@ -254,9 +254,9 @@ def test_repr(self, diffusion_binding): assert 'model=' in repr_str assert 'targets=' in repr_str - -class TestFitBindingWorkflows: - """End-to-end regression tests for the standard ParameterAnalysis workflows.""" + ############# + # Workflows: end-to-end regression tests for the standard ParameterAnalysis workflows + ############# def test_polynomial_targets_gaussian_area(self): # WHEN: fitting a Polynomial to a 'Gaussian area' dataset key diff --git a/tests/unit/easydynamics/analysis/test_parameter_analysis.py b/tests/unit/easydynamics/analysis/test_parameter_analysis.py index 1b6b80178..647c73c22 100644 --- a/tests/unit/easydynamics/analysis/test_parameter_analysis.py +++ b/tests/unit/easydynamics/analysis/test_parameter_analysis.py @@ -9,6 +9,7 @@ import pytest import scipp as sc from easyscience.fitting.multi_fitter import MultiFitter +from easyscience.variable import Parameter import easydynamics as edyn import easydynamics.sample_model as sm @@ -20,6 +21,10 @@ from easydynamics.sample_model.diffusion_model.brownian_translational_diffusion import ( BrownianTranslationalDiffusion, ) +from easydynamics.sample_model.diffusion_model.delta_lorentz import DeltaLorentz +from easydynamics.sample_model.diffusion_model.jump_translational_diffusion import ( + JumpTranslationalDiffusion, +) from easydynamics.utils.fit_target import FitTarget Q = np.array([0.5, 0.8, 1.1, 1.4, 1.7, 2.0]) @@ -574,6 +579,38 @@ def test_plot_no_bindings( # 6. Return value propagated assert result == mock_plot.return_value + def test_plot_with_empty_names_raises_a_clear_error(self, parameter_analysis): + # WHEN / THEN / EXPECT: an empty list is an error, not a bare IndexError + with ( + patch( + 'easydynamics.analysis.parameter_analysis._in_notebook', + return_value=True, + ), + pytest.raises(ValueError, match='names must not be an empty list'), + ): + parameter_analysis.plot(names=[]) + + def test_plot_evaluates_only_the_bindings_being_plotted( + self, parameter_analysis, mock_model_dataset + ): + # WHEN only the first binding's target is requested + parameter_analysis.calculate_model_dataset = MagicMock(return_value=mock_model_dataset) + + # THEN + with ( + patch( + 'easydynamics.analysis.parameter_analysis._in_notebook', + return_value=True, + ), + patch('easydynamics.analysis.parameter_analysis.pp.plot'), + ): + parameter_analysis.plot(names=['parameter1']) + + # EXPECT the diffusion binding is not evaluated for a plot that does not show it + parameter_analysis.calculate_model_dataset.assert_called_once_with([ + parameter_analysis.bindings[0] + ]) + @pytest.mark.parametrize( 'set_pars_none, bindings, expected_exception, match', [ @@ -1139,6 +1176,26 @@ def test_get_xyweight_from_dataset_no_variances(self, parameter_analysis): np.testing.assert_allclose(y, [1.0, 2.0]) np.testing.assert_allclose(w, [1.0, 1.0]) + def test_get_xyweight_from_dataset_no_variances_filters_nan_values(self, parameter_analysis): + # WHEN a dataset without variances contains a NaN value + Q = sc.array(dims=['Q'], values=[0.1, 0.2], unit='1/angstrom') + parameter_analysis.parameters = sc.Dataset( + data={ + 'parameter1': sc.DataArray( + data=sc.array(dims=['Q'], values=[1.0, np.nan], unit='meV'), + coords={'Q': Q}, + ), + } + ) + + # THEN + x, y, w = parameter_analysis._get_xyweight_from_dataset('parameter1') + + # EXPECT the NaN row is filtered like on the with-variances path + np.testing.assert_allclose(x, [0.1]) + np.testing.assert_allclose(y, [1.0]) + np.testing.assert_allclose(w, [1.0]) + def test_get_xyweight_from_dataset_all_nan_variances_raises(self, parameter_analysis): # WHEN Q = sc.array(dims=['Q'], values=[0.1, 0.2], unit='1/angstrom') @@ -1279,6 +1336,134 @@ def test_shrinking_the_targets_also_rebuilds(self): # EXPECT assert len(analysis.fit()) == 1 + def test_swapping_targets_of_the_same_model_rebuilds_the_fitter(self): + # WHEN a binding's targets are swapped in place without changing how many there are: + # the model list stays identical, so a model-only signature would miss the change and + # fit the stale width function against the area data + binding = edyn.FitBinding( + model=sm.BrownianTranslationalDiffusion( + name='Brownian', + lorentzian_name='Lorentzian', + diffusion_coefficient=2.4e-9, + scale=0.5, + ), + targets=['width'], + ) + analysis = edyn.ParameterAnalysis(parameters=make_dataset(), bindings=[binding]) + analysis.fit() + original = analysis._fitter + + # THEN + binding.targets = ['area'] + analysis.fit() + + # EXPECT the fitter was rebuilt for the new target + assert analysis._fitter is not original + + def test_append_binding_invalidates_the_fitter(self, analysis): + # WHEN + original = analysis.fitter + new_binding = edyn.FitBinding( + model=sm.Polynomial(coefficients=[1.0], x_unit='1/angstrom', y_unit='meV'), + targets='Lorentzian width', + ) + + # THEN + analysis.append_binding(new_binding) + + # EXPECT + assert analysis.fitter is not original + + def test_clear_bindings_invalidates_the_fitter(self, analysis): + # WHEN + _ = analysis.fitter + assert analysis._fitter_is_dirty is False + + # THEN + analysis.clear_bindings() + + # EXPECT + assert analysis._fitter_is_dirty is True + + def test_bindings_list_is_copied_from_the_caller(self): + # WHEN a caller passes a list and mutates it afterwards + binding = edyn.FitBinding( + model=sm.Polynomial(coefficients=[1.0], x_unit='1/angstrom', y_unit='meV'), + targets='Lorentzian width', + ) + caller_list = [binding] + analysis = edyn.ParameterAnalysis(parameters=make_dataset(), bindings=caller_list) + + # THEN + caller_list.clear() + + # EXPECT the analysis still holds the binding it was given + assert analysis.bindings == [binding] + + ############# + # The bayesian sampler (the ParameterAnalysis side of the contract) + ############# + + def test_bayesian_returns_the_cached_sampler(self, analysis): + # THEN + sampler = analysis.bayesian + + # EXPECT the same object on second access + assert sampler is analysis.bayesian + + def test_bayesian_is_invalidated_when_the_parameters_change(self, analysis): + # WHEN + sampler = analysis.bayesian + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis.parameters = make_dataset() + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_when_the_bindings_change(self, analysis): + # WHEN + sampler = analysis.bayesian + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis.bindings = analysis.bindings[:1] + + # EXPECT + mock_invalidate.assert_called_once() + + def test_bayesian_is_invalidated_when_a_binding_is_appended(self, analysis): + # WHEN + sampler = analysis.bayesian + new_binding = edyn.FitBinding( + model=sm.Polynomial(coefficients=[1.0], x_unit='1/angstrom', y_unit='meV'), + targets='Lorentzian width', + ) + + # THEN + with patch.object(sampler, 'invalidate') as mock_invalidate: + analysis.append_binding(new_binding) + + # EXPECT + mock_invalidate.assert_called_once() + + def test_missing_parameters_dataset_raises(self): + # WHEN + parameter_analysis = edyn.ParameterAnalysis() + + # THEN EXPECT + with pytest.raises(ValueError, match='No parameters Dataset'): + parameter_analysis.bayesian.sample(samples=10) + + def test_missing_bindings_raises(self): + # WHEN + parameter_analysis = edyn.ParameterAnalysis(parameters=make_dataset()) + + # THEN EXPECT + with pytest.raises(ValueError, match='No fit bindings'): + parameter_analysis.bayesian.sample(samples=10) + ############# # Chain parameters and labels ############# @@ -1345,8 +1530,6 @@ def test_single_binding_keeps_plain_names(self): def test_parameter_from_outside_the_analysis_keeps_its_name(self, analysis): # WHEN a parameter belongs to none of the binding models - from easyscience.variable import Parameter - stranger = Parameter(name='Width line_c0', value=1.0) # THEN EXPECT it is returned unqualified rather than mislabelled @@ -1403,8 +1586,6 @@ def test_colliding_names_with_distinct_models_use_the_display_name(self): def test_ambiguous_name_owned_by_no_model_keeps_its_name(self): # WHEN a parameter shares an ambiguous name but belongs to none of the models - from easyscience.variable import Parameter - analysis = edyn.ParameterAnalysis( parameters=make_dataset(), bindings=[ @@ -1454,8 +1635,6 @@ def _dataset_from_targets(model, Q, unit_overrides=None): def test_delta_lorentz_three_target_simultaneous_fit(self): # WHEN: synthetic width, area, and delta area curves from a known DeltaLorentz - from easydynamics.sample_model.diffusion_model.delta_lorentz import DeltaLorentz - Q = np.linspace(0.4, 2.0, 9) truth = DeltaLorentz(scale=2.0, mean_u_squared=0.3, A_0=0.6, lorentzian_width=0.12) dataset = self._dataset_from_targets(truth, Q) @@ -1473,10 +1652,6 @@ def test_delta_lorentz_three_target_simultaneous_fit(self): def test_jump_diffusion_width_only_fit(self): # WHEN: synthetic widths from a known jump diffusion model - from easydynamics.sample_model.diffusion_model.jump_translational_diffusion import ( - JumpTranslationalDiffusion, - ) - Q = np.linspace(0.4, 2.0, 9) truth = JumpTranslationalDiffusion(diffusion_coefficient=2.4e-9, relaxation_time=2.0) dataset = self._dataset_from_targets(truth, Q) diff --git a/tests/unit/easydynamics/analysis/test_posterior.py b/tests/unit/easydynamics/analysis/test_posterior.py index 436283a4c..40a6a3b66 100644 --- a/tests/unit/easydynamics/analysis/test_posterior.py +++ b/tests/unit/easydynamics/analysis/test_posterior.py @@ -1,6 +1,8 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import warnings + import numpy as np import pytest from easyscience.variable import Parameter @@ -132,13 +134,15 @@ def test_non_numeric_setting_raises(self): suggest_bounds_for_parameters([make_parameter()], n_sigma='wide') -class TestBoundsSuggestionsApply: +class TestBoundsSuggestions: + ############# + # Applying suggestions + ############# + def test_apply_sets_bounds_and_reports_changes(self): - # WHEN + # WHEN nothing has changed yet, since apply has not been called parameter = make_parameter(value=10.0, error=0.5) suggestions = suggest_bounds_for_parameters([parameter]) - - # WHEN nothing has changed until apply is called assert parameter.max == np.inf # THEN @@ -161,6 +165,45 @@ def test_apply_skips_parameters_needing_attention(self): assert changed == [] assert parameter.min == -np.inf + ############# + # Absurd-bounds warning + ############# + + def test_applying_a_wildly_wide_bound_warns(self): + # WHEN a fit returns an enormous uncertainty, which is what a degenerate parameter looks + # like coming out of least squares + parameter = make_parameter(name='Delta area', value=1.0, error=1e9) + suggestions = suggest_bounds_for_parameters([parameter]) + + # THEN EXPECT it is still applied, since it is what the fit implied, but not silently + with pytest.warns(UserWarning, match='far wider than the parameter'): + changed = suggestions.apply() + assert changed == [parameter] + + def test_a_sane_bound_applies_without_warning(self): + # WHEN + parameter = make_parameter(name='sane', value=10.0, error=0.5) + suggestions = suggest_bounds_for_parameters([parameter]) + + # THEN EXPECT + with warnings.catch_warnings(): + warnings.simplefilter('error') + suggestions.apply() + + def test_a_zero_valued_parameter_is_not_called_absurd(self): + # WHEN there is no magnitude to compare the width against + parameter = make_parameter(name='zero', value=0.0, error=1.0) + suggestions = suggest_bounds_for_parameters([parameter]) + + # THEN EXPECT no warning, since the ratio is meaningless rather than alarming + with warnings.catch_warnings(): + warnings.simplefilter('error') + suggestions.apply() + + ############# + # Repr and iteration + ############# + def test_repr_lists_parameters_and_flags_attention(self): # WHEN good = make_parameter(name='good', value=10.0, error=0.5) @@ -347,6 +390,8 @@ def test_unmatched_column_falls_back_to_the_supplied_name(self): assert entry.unit == '' assert np.isnan(entry.value) + +class TestPosteriorSummary: def test_lookup_of_missing_name_raises(self): # WHEN summary = summarize_draws(np.zeros((5, 1)), ['x'], [None]) @@ -366,8 +411,6 @@ def test_repr_contains_the_parameter_name(self): assert 'Gaussian area' in text assert 'median' in text - -class TestPosteriorSummaryContainer: def test_len_and_iteration(self): # WHEN parameters = [make_parameter(name='a'), make_parameter(name='b')] @@ -381,40 +424,3 @@ def test_len_and_iteration(self): def test_repr_with_no_entries(self): # WHEN THEN EXPECT assert 'no parameters' in repr(summarize_draws(np.zeros((3, 0)), [], [])) - - -class TestAbsurdBoundsWarning: - def test_applying_a_wildly_wide_bound_warns(self): - # WHEN a fit returns an enormous uncertainty, which is what a degenerate parameter looks - # like coming out of least squares - parameter = make_parameter(name='Delta area', value=1.0, error=1e9) - suggestions = suggest_bounds_for_parameters([parameter]) - - # THEN EXPECT it is still applied, since it is what the fit implied, but not silently - with pytest.warns(UserWarning, match='far wider than the parameter'): - changed = suggestions.apply() - assert changed == [parameter] - - def test_a_sane_bound_applies_without_warning(self): - # WHEN - parameter = make_parameter(name='sane', value=10.0, error=0.5) - suggestions = suggest_bounds_for_parameters([parameter]) - - # THEN EXPECT - import warnings as warnings_module - - with warnings_module.catch_warnings(): - warnings_module.simplefilter('error') - suggestions.apply() - - def test_a_zero_valued_parameter_is_not_called_absurd(self): - # WHEN there is no magnitude to compare the width against - parameter = make_parameter(name='zero', value=0.0, error=1.0) - suggestions = suggest_bounds_for_parameters([parameter]) - - # THEN EXPECT no warning, since the ratio is meaningless rather than alarming - import warnings as warnings_module - - with warnings_module.catch_warnings(): - warnings_module.simplefilter('error') - suggestions.apply() diff --git a/tests/unit/easydynamics/analysis/test_posterior_labels.py b/tests/unit/easydynamics/analysis/test_posterior_labels.py index 618c8284c..0751c7895 100644 --- a/tests/unit/easydynamics/analysis/test_posterior_labels.py +++ b/tests/unit/easydynamics/analysis/test_posterior_labels.py @@ -11,7 +11,25 @@ def make_parameter(name, unit='meV'): return Parameter(name=name, value=1.0, unit=unit) -class TestLabelling: +class _CountingParameter: + """A Parameter stand-in that counts how often its name is read.""" + + def __init__(self, name, unique_name): + self._name = name + self.unique_name = unique_name + self.name_accesses = 0 + + @property + def name(self): + self.name_accesses += 1 + return self._name + + +class TestParameterLabels: + ############# + # Labelling + ############# + def test_unique_names_are_left_alone(self): # WHEN nothing is ambiguous, a qualifier would only cost width parameters = [make_parameter('area'), make_parameter('width')] @@ -54,8 +72,10 @@ def test_without_a_qualifier_names_stay_bare(self): # EXPECT assert labels.label(first) == 'width' + ############# + # Chain columns + ############# -class TestChainColumns: def test_columns_resolve_by_unique_name(self): # WHEN parameters = [make_parameter('area'), make_parameter('width')] @@ -131,19 +151,25 @@ def test_name_map_records_labels_against_unique_names(self): second.unique_name: 'width (Q_index=1)', } + ############# + # Cost + ############# -class TestCost: def test_labelling_does_not_rescan_per_parameter(self): # WHEN there are many parameters. Computing the name counts per parameter is quadratic, - # which was seconds of work for an analysis with many Q values. - parameters = [make_parameter(f'p{i // 2}') for i in range(400)] + # which was seconds of work for an analysis with many Q values, so labelling one parameter + # must not read every parameter's name again. + parameters = [_CountingParameter(f'p{i // 2}', f'Parameter_{i}') for i in range(400)] labels = ParameterLabels(parameters, qualify=lambda _p: 'q') + for parameter in parameters: + parameter.name_accesses = 0 - # THEN EXPECT labelling all of them stays cheap - import time - - start = time.perf_counter() + # THEN names = [labels.label(p) for p in parameters] - assert time.perf_counter() - start < 0.5 + + # EXPECT a bounded number of name reads per label() call: a quadratic implementation + # recounting the names inside label() would read all 400 names on every call + total_accesses = sum(p.name_accesses for p in parameters) + assert total_accesses <= 4 * len(parameters) assert len(names) == len(parameters) assert np.all([n.endswith('(q)') for n in names]) diff --git a/tests/unit/easydynamics/analysis/test_posterior_sampling.py b/tests/unit/easydynamics/analysis/test_posterior_sampling.py index ab3759f95..e46edc22f 100644 --- a/tests/unit/easydynamics/analysis/test_posterior_sampling.py +++ b/tests/unit/easydynamics/analysis/test_posterior_sampling.py @@ -8,7 +8,9 @@ for PosteriorSampler, and an Analysis for the multi-Q subclass. """ +import json import types +import warnings from types import SimpleNamespace from unittest.mock import MagicMock from unittest.mock import patch @@ -18,11 +20,12 @@ import pytest import scipp as sc from easyscience.fitting import AvailableMinimizers -from easyscience.fitting.multi_fitter import MultiFitter from easyscience.variable import Parameter mpl.use('Agg') +import matplotlib.pyplot as plt + import easydynamics as edyn import easydynamics.sample_model as sm from easydynamics.analysis.analysis1d import Analysis1d @@ -560,8 +563,6 @@ def test_save_without_a_chain_raises(self, analysis): def test_save_writes_the_parameter_name_sidecar(self, analysis, tmp_path): # WHEN - import json - bound_all(analysis) # THEN @@ -798,10 +799,6 @@ def test_load_chain_uses_the_sidecar_when_present(self, analysis, tmp_path): def test_trace_and_corner_render_from_a_chain(self, analysis): # WHEN - import matplotlib as mpl - import matplotlib.pyplot as plt - - mpl.use('Agg') bound_all(analysis) n_parameters = len(analysis.get_free_parameters()) @@ -1219,22 +1216,6 @@ def mutate(**_kwargs): # EXPECT assert [float(p.value) for p in parameters] == pytest.approx(before) - def test_missing_parameters_dataset_raises(self): - # WHEN - parameter_analysis = edyn.ParameterAnalysis() - - # THEN EXPECT - with pytest.raises(ValueError, match='No parameters Dataset'): - parameter_analysis.bayesian.sample(samples=10) - - def test_missing_bindings_raises(self): - # WHEN - parameter_analysis = edyn.ParameterAnalysis(parameters=make_dataset()) - - # THEN EXPECT - with pytest.raises(ValueError, match='No fit bindings'): - parameter_analysis.bayesian.sample(samples=10) - class TestMultiQPosteriorSampler: ############# @@ -1325,13 +1306,6 @@ def test_refreshes_every_convolver_before_sampling(self, multi_q_analysis): # EXPECT the sampler sees the same prepared convolvers a simultaneous fit would assert all(not a._convolver_is_dirty for a in multi_q_analysis.analysis_list) - def test_uses_a_multifitter(self, multi_q_analysis): - # WHEN - - # EXPECT - assert isinstance(multi_q_analysis.fitter, MultiFitter) - assert len(multi_q_analysis.fitter.fit_object) == len(Q_VALUES) - ############# # Independent sampling ############# @@ -1995,8 +1969,6 @@ class warnings_as_errors: """Context manager asserting that no UserWarning is emitted inside the block.""" def __enter__(self): - import warnings - self._ctx = warnings.catch_warnings(record=True) self._caught = self._ctx.__enter__() warnings.simplefilter('always') diff --git a/tests/unit/easydynamics/base_classes/test_easydynamics_list.py b/tests/unit/easydynamics/base_classes/test_easydynamics_list.py index 3d81ec60e..aa7e4d37f 100644 --- a/tests/unit/easydynamics/base_classes/test_easydynamics_list.py +++ b/tests/unit/easydynamics/base_classes/test_easydynamics_list.py @@ -257,3 +257,118 @@ def test_getitem_invalid_type(self, easy_dynamics_list): # WHEN THEN EXPECT with pytest.raises(TypeError, match=r'Index must be an int, slice, or str'): easy_dynamics_list[1.5] + + ############# + # Item assignment + ############# + + def test_setitem(self, easy_dynamics_list): + """Test assigning an item by index replaces it and nothing else.""" + # WHEN + new_gaussian = Gaussian(name='ReplacementGaussian') + + # THEN + easy_dynamics_list[0] = new_gaussian + + # EXPECT + assert easy_dynamics_list[0] is new_gaussian + assert len(easy_dynamics_list) == 2 + + def test_setitem_invalid_type_raises(self, easy_dynamics_list): + # WHEN THEN EXPECT + with pytest.raises(TypeError): + easy_dynamics_list[0] = 'Not a ModelComponent' + + def test_setitem_repeated_component_warns(self, easy_dynamics_list): + """Test that item assignment warns and ignores like append/insert do.""" + # WHEN THEN EXPECT assigning an item already in the list warns and is ignored + with pytest.warns(UserWarning, match=r'already in EasyDynamicsList'): + easy_dynamics_list[1] = easy_dynamics_list[0] + + assert easy_dynamics_list[1] is not easy_dynamics_list[0] + + ############# + # Versioning + ############# + + def test_version_starts_at_zero(self, easy_dynamics_list): + # WHEN a freshly constructed list, even with initial items + # THEN EXPECT version is 0 + assert easy_dynamics_list.version == 0 + + def test_version_is_read_only(self, easy_dynamics_list): + # WHEN THEN EXPECT + with pytest.raises(AttributeError): + easy_dynamics_list.version = 5 + + def test_version_bumps_on_every_mutator(self, easy_dynamics_list): + """Every mutating operation increments version; reads do not.""" + # WHEN + version = easy_dynamics_list.version + + # THEN append + easy_dynamics_list.append(Gaussian(name='V1')) + # EXPECT + assert easy_dynamics_list.version == version + 1 + + # THEN insert + easy_dynamics_list.insert(0, Gaussian(name='V2')) + # EXPECT + assert easy_dynamics_list.version == version + 2 + + # THEN extend (one bump per item) + easy_dynamics_list.extend([Gaussian(name='V3'), Gaussian(name='V4')]) + # EXPECT + assert easy_dynamics_list.version == version + 4 + + # THEN item assignment + easy_dynamics_list[0] = Gaussian(name='V5') + # EXPECT + assert easy_dynamics_list.version == version + 5 + + # THEN pop by index and by name + easy_dynamics_list.pop(0) + easy_dynamics_list.pop('V1') + # EXPECT + assert easy_dynamics_list.version == version + 7 + + # THEN remove and del + item = easy_dynamics_list[0] + easy_dynamics_list.remove(item) + del easy_dynamics_list[0] + # EXPECT + assert easy_dynamics_list.version == version + 9 + + # THEN sort + easy_dynamics_list.sort(key=lambda c: c.name) + # EXPECT + assert easy_dynamics_list.version == version + 10 + + # THEN clear + n_items = len(easy_dynamics_list) + easy_dynamics_list.clear() + # EXPECT one bump per removed item, and reading version mutates nothing + assert easy_dynamics_list.version == version + 10 + n_items + assert easy_dynamics_list.version == version + 10 + n_items + + def test_version_does_not_bump_on_ignored_duplicate(self, easy_dynamics_list): + # WHEN + version = easy_dynamics_list.version + + # THEN an insert that is ignored because the item is already in the list + with pytest.warns(UserWarning, match=r'already in EasyDynamicsList'): + easy_dynamics_list.insert(1, easy_dynamics_list[0]) + + # EXPECT no mutation happened, so no version bump + assert easy_dynamics_list.version == version + + def test_version_does_not_bump_on_failed_mutation(self, easy_dynamics_list): + # WHEN + version = easy_dynamics_list.version + + # THEN EXPECT failed mutations leave the version unchanged + with pytest.raises(TypeError): + easy_dynamics_list.append('Not a ModelComponent') + with pytest.raises(KeyError): + easy_dynamics_list.pop('Nonexistent') + assert easy_dynamics_list.version == version diff --git a/tests/unit/easydynamics/base_classes/test_name_mixin.py b/tests/unit/easydynamics/base_classes/test_name_mixin.py index 5ea931633..1913e519e 100644 --- a/tests/unit/easydynamics/base_classes/test_name_mixin.py +++ b/tests/unit/easydynamics/base_classes/test_name_mixin.py @@ -2,6 +2,8 @@ # SPDX-License-Identifier: BSD-3-Clause import pytest +from easyscience import global_object +from easyscience.base_classes.new_base import NewBase from easydynamics.base_classes.name_mixin import NameMixin @@ -62,3 +64,19 @@ def test_name_setter_invalid_type(self, name_mixin, invalid_name): # WHEN THEN EXPECT with pytest.raises(TypeError, match=r'Name must be a string.'): name_mixin.name = invalid_name + + def test_invalid_name_fails_before_global_registration(self): + """Regression: name validation must run before the parent registers the object.""" + + # WHEN a class whose MRO reaches the registering NewBase through NameMixin + class _RegisteredWithName(NameMixin, NewBase): + pass + + vertices_before = set(global_object.map.vertices()) + + # THEN EXPECT construction fails on the invalid name + with pytest.raises(TypeError, match=r'Name must be a string'): + _RegisteredWithName(name=123) + + # EXPECT no half-constructed object was registered in the global map + assert set(global_object.map.vertices()) == vertices_before diff --git a/tests/unit/easydynamics/convolution/test_convolution.py b/tests/unit/easydynamics/convolution/test_convolution.py index 15305ae54..ea50c2f0e 100644 --- a/tests/unit/easydynamics/convolution/test_convolution.py +++ b/tests/unit/easydynamics/convolution/test_convolution.py @@ -8,6 +8,7 @@ import numpy as np import pytest import scipp as sc +from easyscience.variable import Parameter from easydynamics.convolution.analytical_convolution import AnalyticalConvolution from easydynamics.convolution.convolution import Convolution @@ -20,6 +21,7 @@ from easydynamics.sample_model import Polynomial from easydynamics.sample_model import Voigt from easydynamics.sample_model.component_collection import ComponentCollection +from easydynamics.settings.detailed_balance_settings import DetailedBalanceSettings class TestConvolution: @@ -402,7 +404,7 @@ def test_check_if_pair_is_analytic(self, default_convolution, function1, functio def test_check_if_pair_is_analytic_raises_with_delta_in_resolution(self, default_convolution): """ - Test that _check_if_pair_is_analytic raises TypeError when + Test that _check_if_pair_is_analytic raises ValueError when resolution component is DeltaFunction. """ # WHEN @@ -412,7 +414,7 @@ def test_check_if_pair_is_analytic_raises_with_delta_in_resolution(self, default # THEN EXPECT with pytest.raises( - TypeError, + ValueError, match='This is not supported', ): conv._check_if_pair_is_analytic( @@ -632,6 +634,251 @@ def test_convert_y_unit_propagates_to_sub_convolvers(self): assert conv.y_unit == '1/eV' assert conv._analytical_convolver._y_unit == '1/eV' + ############# + # Plan invalidation regressions + ############# + + def test_invalidate_plan_on_change_names_are_real_attributes(self, default_convolution): + "Regression: the tracked-attribute set used to contain names that never exist" + # WHEN THEN EXPECT every tracked name is an actual attribute of a built convolver + for name in Convolution._invalidate_plan_on_change: + assert hasattr(default_convolution, name), name + + def test_in_place_sample_append_contributes_to_convolution(self): + "Regression: appending to the live sample collection used to leave output unchanged" + # WHEN a convolver that has already produced output + energy = np.linspace(-10, 10, 1001) + sample_components = ComponentCollection( + components=[Gaussian(name='G', area=2.0, center=0.1, width=0.4)] + ) + resolution_components = ComponentCollection( + components=[Gaussian(name='R', area=1.0, center=0.0, width=0.5)] + ) + conv = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + ) + result_before = conv.convolution() + + # THEN mutating the live sample collection in place + conv.sample_components.append_component( + Lorentzian(name='L', area=1.0, center=0.0, width=0.3) + ) + result_after = conv.convolution() + + # EXPECT the new component contributes, matching a freshly built convolver + fresh = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + ) + assert not np.allclose(result_after, result_before) + np.testing.assert_allclose(result_after, fresh.convolution(), rtol=1e-10) + + def test_in_place_resolution_append_contributes_to_convolution(self): + "Regression: appending to the live resolution collection used to leave output unchanged" + # WHEN a convolver that has already produced output + energy = np.linspace(-10, 10, 1001) + sample_components = ComponentCollection( + components=[Gaussian(name='G', area=2.0, center=0.1, width=0.4)] + ) + resolution_components = ComponentCollection( + components=[Gaussian(name='R', area=1.0, center=0.0, width=0.5)] + ) + conv = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + ) + result_before = conv.convolution() + + # THEN mutating the live resolution collection in place + conv.resolution_components.append_component( + Gaussian(name='R2', area=0.5, center=0.0, width=0.2) + ) + result_after = conv.convolution() + + # EXPECT the new component contributes, matching a freshly built convolver + fresh = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + ) + assert not np.allclose(result_after, result_before) + np.testing.assert_allclose(result_after, fresh.convolution(), rtol=1e-10) + + def test_detailed_balance_toggle_changes_output(self): + "Regression: toggling use_detailed_balance after construction used to be ignored" + # WHEN a convolver built with detailed balance off + energy = np.linspace(-10, 10, 1001) + sample_components = ComponentCollection( + components=[Lorentzian(name='L', area=2.0, center=0.0, width=0.4)] + ) + resolution_components = ComponentCollection( + components=[Gaussian(name='R', area=1.0, center=0.0, width=0.5)] + ) + conv = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + temperature=300.0, + detailed_balance_settings=DetailedBalanceSettings(use_detailed_balance=False), + ) + result_off = conv.convolution() + + # THEN toggling detailed balance on after construction + conv.detailed_balance_settings.use_detailed_balance = True + result_on = conv.convolution() + + # EXPECT the output changes and matches a convolver built with detailed balance on + fresh = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + temperature=300.0, + detailed_balance_settings=DetailedBalanceSettings(use_detailed_balance=True), + ) + assert not np.allclose(result_on, result_off) + np.testing.assert_allclose(result_on, fresh.convolution(), rtol=1e-10) + + def test_energy_offset_rebind_reaches_sub_convolvers(self): + "Regression: rebinding energy_offset used to leave sub-convolvers on the old Parameter" + # WHEN a convolver with analytical, numerical and delta components and offset 0 + energy = np.linspace(-10, 10, 1001) + sample_components = ComponentCollection( + components=[ + Gaussian(name='G', area=2.0, center=0.1, width=0.4), + DampedHarmonicOscillator(name='DHO', area=2.0, center=1.0, width=0.1), + DeltaFunction(name='D', area=1.0, center=0.3), + ] + ) + resolution_components = ComponentCollection( + components=[Gaussian(name='R', area=1.0, center=0.0, width=0.5)] + ) + conv = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + ) + result_before = conv.convolution() + + # THEN rebinding the offset to a brand-new Parameter + conv.energy_offset = Parameter(name='energy_offset', value=1.0, unit='meV') + result_after = conv.convolution() + + # EXPECT every path (analytical, numerical, delta) sees the new offset, matching a + # freshly built convolver + fresh = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + energy_offset=1.0, + ) + assert not np.allclose(result_after, result_before) + np.testing.assert_allclose(result_after, fresh.convolution(), rtol=1e-10) + + ############# + # Dispatch and validation regressions + ############# + + def test_subclass_of_analytical_component_convolves_like_base(self): + "Regression: a Lorentzian subclass was routed analytically but rejected by dispatch" + + # WHEN a subclass of Lorentzian in the sample model + class MyLorentzian(Lorentzian): + pass + + energy = np.linspace(-10, 10, 1001) + conv = Convolution( + energy=energy, + sample_components=MyLorentzian(name='MyL', area=2.0, center=0.1, width=0.4), + resolution_components=Gaussian(name='R', area=1.0, center=0.0, width=0.5), + ) + + # THEN it is routed to the analytical convolver and convolved with the base rules + result = conv.convolution() + + # EXPECT + assert len(conv._analytical_sample_components) == 1 + reference = Convolution( + energy=energy, + sample_components=Lorentzian(name='L', area=2.0, center=0.1, width=0.4), + resolution_components=Gaussian(name='R2', area=1.0, center=0.0, width=0.5), + ) + np.testing.assert_allclose(result, reference.convolution(), rtol=1e-10) + + def test_empty_resolution_raises(self): + "Regression: an empty resolution used to silently produce zeros" + # WHEN THEN EXPECT at construction + with pytest.raises(ValueError, match=r'resolution_components is empty'): + Convolution( + energy=np.linspace(-10, 10, 101), + sample_components=Gaussian(name='G', area=1.0, center=0.0, width=0.4), + resolution_components=ComponentCollection(), + ) + + def test_emptying_resolution_in_place_raises_on_next_convolution(self, default_convolution): + # WHEN the live resolution collection is emptied after construction + conv = default_convolution + conv.resolution_components.pop('GaussianRes') + + # THEN EXPECT the next convolution rebuilds the plan and refuses to silently + # return zeros + with pytest.raises(ValueError, match=r'resolution_components is empty'): + conv.convolution() + + ############# + # Registry and label housekeeping + ############# + + def test_plan_rebuilds_do_not_leak_registry_entries(self, default_convolution): + "Regression: every plan rebuild used to register new objects in the global map forever" + # WHEN a convolver that has built its plan at least once + conv = default_convolution + conv.convolution() + vertices_before = len(conv._global_object.map.vertices()) + + # THEN forcing several full plan rebuilds + for _ in range(3): + conv._plan_seen_version = None + conv.convolution() + + # EXPECT the global map did not grow + assert len(conv._global_object.map.vertices()) == vertices_before + + def test_convert_y_unit_updates_plan_collection_labels(self): + "Regression: plan-collection y_unit labels used to stay stale until the next rebuild" + # WHEN a convolver with analytical, numerical and delta components in 1/meV + energy = np.linspace(-10, 10, 1001) + sample_components = ComponentCollection( + components=[ + Gaussian(name='G', area=1.0, center=0.0, width=0.4, y_unit='1/meV'), + DampedHarmonicOscillator( + name='DHO', area=1.0, center=1.0, width=0.1, y_unit='1/meV' + ), + DeltaFunction(name='D', area=1.0, center=0.0, y_unit='1/meV'), + ], + y_unit='1/meV', + ) + resolution_components = ComponentCollection( + components=[Gaussian(name='R', area=1.0, center=0.0, width=0.5)] + ) + conv = Convolution( + energy=energy, + sample_components=sample_components, + resolution_components=resolution_components, + y_unit='1/meV', + ) + + # THEN + conv.convert_y_unit('1/eV') + + # EXPECT the plan collections' labels follow without waiting for a rebuild + assert conv._analytical_sample_components.y_unit == '1/eV' + assert conv._numerical_sample_components.y_unit == '1/eV' + assert conv._delta_sample_components.y_unit == '1/eV' + def test_convert_y_unit_propagates_to_numerical_convolver(self): # WHEN: a DHO sample component forces a numerical convolver energy = np.linspace(-10, 10, 5001) diff --git a/tests/unit/easydynamics/convolution/test_convolution_base.py b/tests/unit/easydynamics/convolution/test_convolution_base.py index 7bc26fff7..fedf214e0 100644 --- a/tests/unit/easydynamics/convolution/test_convolution_base.py +++ b/tests/unit/easydynamics/convolution/test_convolution_base.py @@ -8,6 +8,7 @@ from scipp import UnitError from easydynamics.convolution.convolution_base import ConvolutionBase +from easydynamics.sample_model import DeltaFunction from easydynamics.sample_model import Gaussian from easydynamics.sample_model.component_collection import ComponentCollection @@ -403,3 +404,85 @@ def test_convert_y_unit_without_sample_components(self): # EXPECT assert cb.y_unit == '1/meV' + + ############# + # Unit-consistency validation + ############# + + def test_energy_setter_scipp_with_matching_unit(self, convolution_base): + # WHEN + new_energy = sc.array(dims=['energy'], values=np.linspace(-3, 3, 7), unit='meV') + + # THEN + convolution_base.energy = new_energy + + # EXPECT: accepted and x_unit stays a str + assert sc.identical(convolution_base.energy, new_energy) + assert isinstance(convolution_base.x_unit, str) + assert convolution_base.x_unit == 'meV' + + def test_energy_setter_scipp_with_different_unit_raises(self, convolution_base): + "Regression: a mismatched scipp energy used to silently overwrite x_unit with sc.Unit" + # WHEN + new_energy = sc.array(dims=['energy'], values=np.linspace(-3, 3, 7), unit='ueV') + + # THEN EXPECT: unit changes must go through convert_x_unit + with pytest.raises(ValueError, match=r'Use convert_x_unit'): + convolution_base.energy = new_energy + + # EXPECT: nothing changed + assert convolution_base.x_unit == 'meV' + assert np.allclose(convolution_base.energy.values, np.linspace(-10, 10, 100)) + + def test_init_sample_components_x_unit_mismatch_raises(self): + # WHEN sample components in ueV but the convolver in meV + sample = ComponentCollection(components=Gaussian(name='G', x_unit='ueV'), x_unit='ueV') + + # THEN EXPECT + with pytest.raises(ValueError, match=r'sample_components has x_unit'): + ConvolutionBase( + energy=np.linspace(-10, 10, 100), + sample_components=sample, + resolution_components=ComponentCollection(), + x_unit='meV', + ) + + def test_init_resolution_components_x_unit_mismatch_raises(self): + # WHEN resolution components in ueV but the convolver in meV + resolution = ComponentCollection(components=Gaussian(name='R', x_unit='ueV'), x_unit='ueV') + + # THEN EXPECT + with pytest.raises(ValueError, match=r'resolution_components has x_unit'): + ConvolutionBase( + energy=np.linspace(-10, 10, 100), + sample_components=ComponentCollection(), + resolution_components=resolution, + x_unit='meV', + ) + + ############# + # Delta functions in the resolution + ############# + + def test_init_with_delta_in_resolution_raises(self): + # WHEN + resolution = ComponentCollection(components=DeltaFunction(name='D')) + + # THEN EXPECT + with pytest.raises(ValueError, match=r'delta functions'): + ConvolutionBase( + energy=np.linspace(-10, 10, 100), + sample_components=ComponentCollection(), + resolution_components=resolution, + ) + + def test_resolution_components_setter_with_delta_raises(self, convolution_base): + # WHEN + resolution = ComponentCollection(components=DeltaFunction(name='D')) + + # THEN EXPECT + with pytest.raises(ValueError, match=r'delta functions'): + convolution_base.resolution_components = resolution + + # EXPECT: the previous resolution model is kept + assert convolution_base.resolution_components is not resolution diff --git a/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py b/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py index 8fd4722d9..3d80f983b 100644 --- a/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py +++ b/tests/unit/easydynamics/convolution/test_numerical_convolution_base.py @@ -9,6 +9,7 @@ from easydynamics.convolution.energy_grid import EnergyGrid from easydynamics.convolution.numerical_convolution_base import NumericalConvolutionBase from easydynamics.sample_model import Gaussian +from easydynamics.sample_model import Voigt from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.settings.convolution_settings import ConvolutionSettings from easydynamics.settings.detailed_balance_settings import DetailedBalanceSettings @@ -621,6 +622,149 @@ def test_check_width_no_warnings(self, default_numerical_convolution_base): model_name='ComponentCollection', ) + def test_init_with_none_sample_components_raises(self): + # WHEN THEN EXPECT: None components fail early with a clear error + with pytest.raises(TypeError, match=r'sample_components must be .* not None'): + NumericalConvolutionBase( + energy=np.linspace(-10, 10, 101), + sample_components=None, + resolution_components=ComponentCollection(display_name='ResolutionModel'), + ) + + def test_init_with_none_resolution_components_raises(self): + # WHEN THEN EXPECT: None components fail early with a clear error + with pytest.raises(TypeError, match=r'resolution_components must be .* not None'): + NumericalConvolutionBase( + energy=np.linspace(-10, 10, 101), + sample_components=ComponentCollection(display_name='ComponentCollection'), + resolution_components=None, + ) + + @pytest.mark.parametrize('upsample_factor', [None, 5], ids=['no_upsampling', 'upsample_5']) + def test_single_point_energy_raises_clear_error(self, upsample_factor): + """ + Regression: a single energy point used to hit an IndexError (upsample None) or + silently return zeros (default path) instead of the intended ValueError. + """ + # WHEN THEN EXPECT (the grid is built eagerly during construction) + with pytest.raises(ValueError, match=r'at least two points'): + NumericalConvolutionBase( + energy=np.array([1.0]), + sample_components=ComponentCollection(display_name='ComponentCollection'), + resolution_components=ComponentCollection(display_name='ResolutionModel'), + convolution_settings=ConvolutionSettings(upsample_factor=upsample_factor), + ) + + def test_extension_factor_setter_accepts_none(self, default_numerical_convolution_base): + # WHEN + default_numerical_convolution_base.upsample_factor = None + + # THEN + default_numerical_convolution_base.extension_factor = None + + # EXPECT + assert default_numerical_convolution_base.extension_factor is None + + def test_check_width_thresholds_covers_voigt_widths(self, default_numerical_convolution_base): + """ + Regression: width warnings used to gate on 'width' only, silently skipping Voigt + components with gaussian_width/lorentzian_width. + """ + # WHEN a Voigt with one very narrow and one very wide width + voigt = Voigt( + name='NarrowWideVoigt', + area=1.0, + center=0.0, + gaussian_width=1e-6, + lorentzian_width=15.0, + ) + + # THEN EXPECT both widths trigger their warning + with pytest.warns(UserWarning) as record: + default_numerical_convolution_base._check_width_thresholds( + model=voigt, + model_name='sample model', + ) + messages = [str(w.message) for w in record] + assert any('gaussian width' in m and 'upsample_factor' in m for m in messages) + assert any('lorentzian width' in m and 'extension_factor' in m for m in messages) + + ############# + # Plan invalidation + ############# + + def test_detailed_balance_flag_toggle_invalidates_plan( + self, default_numerical_convolution_base + ): + "Regression: toggling detailed balance flags used to be silently ignored" + # WHEN a convolver with a current plan + conv = default_numerical_convolution_base + conv._mark_convolution_plan_current() + assert conv._convolution_plan_is_current() is True + + # THEN + conv.detailed_balance_settings.use_detailed_balance = False + + # EXPECT + assert conv._convolution_plan_is_current() is False + + def test_detailed_balance_settings_rebind_invalidates_plan( + self, default_numerical_convolution_base + ): + # WHEN a convolver with a current plan + conv = default_numerical_convolution_base + conv._mark_convolution_plan_current() + assert conv._convolution_plan_is_current() is True + + # THEN + conv.detailed_balance_settings = DetailedBalanceSettings() + + # EXPECT + assert conv._convolution_plan_is_current() is False + + def test_in_place_collection_mutation_invalidates_plan( + self, default_numerical_convolution_base + ): + "Regression: appending to a live collection used to leave the plan current" + # WHEN a convolver with a current plan + conv = default_numerical_convolution_base + conv._mark_convolution_plan_current() + assert conv._convolution_plan_is_current() is True + + # THEN + conv.sample_components.append_component(Gaussian(name='LiveGaussian')) + + # EXPECT + assert conv._convolution_plan_is_current() is False + + def test_energy_offset_rebind_invalidates_plan(self, default_numerical_convolution_base): + "Regression: rebinding energy_offset to a new Parameter used to cause split-brain" + # WHEN a convolver with a current plan + conv = default_numerical_convolution_base + conv._mark_convolution_plan_current() + + # THEN a numeric assignment mutates the shared Parameter: plan stays current + conv.energy_offset = 1.5 + assert conv._convolution_plan_is_current() is True + + # THEN rebinding to a new Parameter object invalidates the plan + conv.energy_offset = Parameter(name='energy_offset', value=1.5, unit='meV') + + # EXPECT + assert conv._convolution_plan_is_current() is False + + def test_convert_x_unit_invalidates_plan(self, default_numerical_convolution_base): + # WHEN a convolver with a current plan + conv = default_numerical_convolution_base + conv._mark_convolution_plan_current() + + # THEN + conv.convert_x_unit('eV') + + # EXPECT + assert conv._convolution_plan_is_current() is False + assert conv.x_unit == 'eV' + def test_repr(self, default_numerical_convolution_base): """ Test the __repr__ method of NumericalConvolutionBase. @@ -646,14 +790,14 @@ def test_repr(self, default_numerical_convolution_base): assert 'temperature=None' in repr_str assert 'normalize_detailed_balance=True' in repr_str - -def test_create_energy_grid_raises_when_extension_factor_none_with_upsampling(): - # GIVEN upsampling enabled but no extension_factor, the dense energy grid cannot be built - # WHEN THEN EXPECT (the grid is built eagerly during construction) - with pytest.raises(ValueError, match=r'extension_factor must be a number'): - NumericalConvolutionBase( - energy=np.linspace(-10, 10, 101), - sample_components=ComponentCollection(display_name='ComponentCollection'), - resolution_components=ComponentCollection(display_name='ResolutionModel'), - convolution_settings=ConvolutionSettings(upsample_factor=5, extension_factor=None), - ) + def test_create_energy_grid_raises_when_extension_factor_none_with_upsampling(self): + # WHEN upsampling is enabled but there is no extension_factor, the dense energy grid + # cannot be built + # THEN EXPECT (the grid is built eagerly during construction) + with pytest.raises(ValueError, match=r'extension_factor must be a number'): + NumericalConvolutionBase( + energy=np.linspace(-10, 10, 101), + sample_components=ComponentCollection(display_name='ComponentCollection'), + resolution_components=ComponentCollection(display_name='ResolutionModel'), + convolution_settings=ConvolutionSettings(upsample_factor=5, extension_factor=None), + ) diff --git a/tests/unit/easydynamics/experiment/test_experiment.py b/tests/unit/easydynamics/experiment/test_experiment.py index 1aebc0436..cae367c77 100644 --- a/tests/unit/easydynamics/experiment/test_experiment.py +++ b/tests/unit/easydynamics/experiment/test_experiment.py @@ -255,6 +255,19 @@ def test_rebin_with_bin_edge_coordinate(self): assert rebinned_data.sizes['Q'] == 10 assert rebinned_data.sizes['energy'] == 7 + def test_rebin_does_not_mutate_the_callers_dimensions_dict(self, experiment): + "Regression: rebin must not write int-converted values back into the caller's dict" + # WHEN + dimensions = {'Q': 6.0, 'energy': 7} + original = dict(dimensions) + + # THEN + experiment.rebin(dimensions) + + # EXPECT the caller's dict is unchanged (6.0 not silently replaced by 6) + assert dimensions == original + assert isinstance(dimensions['Q'], float) + def test_rebin_no_data_raises(self): "Test rebinning data when no data is present" # WHEN diff --git a/tests/unit/easydynamics/sample_model/components/test_damped_harmonic_oscillator.py b/tests/unit/easydynamics/sample_model/components/test_damped_harmonic_oscillator.py index 12f73235f..7525f1585 100644 --- a/tests/unit/easydynamics/sample_model/components/test_damped_harmonic_oscillator.py +++ b/tests/unit/easydynamics/sample_model/components/test_damped_harmonic_oscillator.py @@ -143,6 +143,24 @@ def test_width_must_be_positive(self, dho: DampedHarmonicOscillator): with pytest.raises(ValueError, match='width must be positive'): dho.width = -0.5 + def test_area_setter_out_of_bounds_raises(self, dho: DampedHarmonicOscillator): + # WHEN the fixture's area was created non-negative, so it carries min=0 + original_area = dho.area.value + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + dho.area = -1.0 + assert dho.area.value == pytest.approx(original_area) + + def test_width_setter_below_minimum_raises(self, dho: DampedHarmonicOscillator): + # WHEN the width parameter carries an absolute minimum (1e-10) + original_width = dho.width.value + + # THEN EXPECT a tiny positive width below the bound raises instead of being clamped + with pytest.raises(ValueError, match='violates the parameter bounds'): + dho.width = 1e-12 + assert dho.width.value == pytest.approx(original_width) + def test_evaluate(self, dho: DampedHarmonicOscillator): # WHEN x = np.array([0.0, 1.5, 3.0]) diff --git a/tests/unit/easydynamics/sample_model/components/test_delta_function.py b/tests/unit/easydynamics/sample_model/components/test_delta_function.py index 85fbf1f6f..7972812fd 100644 --- a/tests/unit/easydynamics/sample_model/components/test_delta_function.py +++ b/tests/unit/easydynamics/sample_model/components/test_delta_function.py @@ -136,6 +136,23 @@ def test_evaluate_unsorted_grid(self): # EXPECT: spike at x=0 with bin width from the sorted grid [0, 1, 2] -> 1.0 np.testing.assert_allclose(result, [1.0, 0.0, 0.0]) + def test_evaluate_single_point_raises(self): + # WHEN: a single x value defines no bin width for the area / bin_width spike + delta = DeltaFunction(area=1.0) + + # THEN EXPECT + with pytest.raises(ValueError, match='single x value'): + delta.evaluate(0.0) + + def test_area_setter_out_of_bounds_raises(self, delta_function: DeltaFunction): + # WHEN the fixture's area was created non-negative, so it carries min=0 + original_area = delta_function.area.value + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + delta_function.area = -1.0 + assert delta_function.area.value == pytest.approx(original_area) + def test_evaluate_out_of_bounds(self, delta_function: DeltaFunction): # WHEN x = np.linspace(1, 2, 100) # center is at 0.5, so out of bounds diff --git a/tests/unit/easydynamics/sample_model/components/test_expression_component.py b/tests/unit/easydynamics/sample_model/components/test_expression_component.py index ba34ff5ff..4243ca8dd 100644 --- a/tests/unit/easydynamics/sample_model/components/test_expression_component.py +++ b/tests/unit/easydynamics/sample_model/components/test_expression_component.py @@ -14,6 +14,9 @@ from easydynamics.sample_model import Gaussian from easydynamics.sample_model import Lorentzian +GAUSSIAN_EXPRESSION = 'A / (sigma*sqrt(2*pi)) * exp(-(x - x0)**2 / (2*sigma**2))' +GAUSSIAN_UNITS = {'A': 'meV', 'x0': 'meV', 'sigma': 'meV'} + class TestExpressionComponent: @pytest.fixture @@ -113,6 +116,15 @@ def test_invalid_function_raises(self): with pytest.raises(ValueError, match='Unsupported function'): ExpressionComponent('A * unknown_func(x)') + @pytest.mark.parametrize('colliding', ['name', 'expression', 'x_unit']) + def test_symbol_colliding_with_attribute_raises(self, colliding): + # WHEN a symbol shadows an existing class attribute, attribute reads would resolve to + # the class attribute while writes hit the parameter, silently diverging + + # THEN EXPECT the collision is rejected at construction + with pytest.raises(ValueError, match='collides with an existing attribute'): + ExpressionComponent(f'{colliding} * x', parameters={colliding: 1.0}) + @pytest.mark.parametrize( 'parameters', [ @@ -359,22 +371,17 @@ def test_erf(self): expected = np.array([-0.84270079, 0.0, 0.84270079]) # erf(-1), erf(0), erf(1) np.testing.assert_allclose(result, expected, rtol=1e-5) + def test_evaluate_raises_when_input_unit_differs_from_x_unit(self): + # WHEN an ExpressionComponent with x_unit meV + expr = ExpressionComponent('A * x', parameters={'A': 2.0}, x_unit='meV') + x = sc.array(dims=['x'], values=[1.0, 2.0], unit='ueV') + # THEN EXPECT a UnitError when evaluating with x in a different unit + with pytest.raises(sc.UnitError, match=r'cannot auto-convert its parameters'): + expr.evaluate(x) -def test_evaluate_raises_when_input_unit_differs_from_x_unit(): - # GIVEN an ExpressionComponent with x_unit meV - expr = ExpressionComponent('A * x', parameters={'A': 2.0}, x_unit='meV') - x = sc.array(dims=['x'], values=[1.0, 2.0], unit='ueV') - # WHEN evaluating with x in a different unit THEN EXPECT a UnitError - with pytest.raises(sc.UnitError, match=r'cannot auto-convert its parameters'): - expr.evaluate(x) - - -GAUSSIAN_EXPRESSION = 'A / (sigma*sqrt(2*pi)) * exp(-(x - x0)**2 / (2*sigma**2))' -GAUSSIAN_UNITS = {'A': 'meV', 'x0': 'meV', 'sigma': 'meV'} - - -class TestExpressionComponentUnitCorrectness: - """Compare unit-aware expressions against the built-in components.""" + ############# + # Unit correctness: comparisons against the built-in components + ############# @pytest.fixture def gaussian_expr(self): @@ -459,8 +466,10 @@ def test_unit_agnostic_expression_does_not_warn(self): x_unit='meV', ) + ############# + # Output unit + ############# -class TestExpressionComponentOutputUnit: def test_output_unit_gaussian_is_dimensionless(self): # WHEN: area in meV divided by sigma in meV expr = ExpressionComponent( @@ -632,12 +641,71 @@ def test_set_unit_warns_when_breaking_consistency(self): x_unit='meV', ) - # THEN EXPECT: relabelling A breaks the output unit + # THEN EXPECT: relabelling A to an incompatible dimension breaks the output unit with pytest.warns(UserWarning, match='does not match'): + expr.set_unit('A', 's/meV') + + def test_set_unit_to_a_convertible_output_unit_rescales_instead_of_warning(self): + # WHEN: a consistent expression whose output stays dimensionless-compatible + expr = ExpressionComponent( + 'A * (x - x0)', + parameters={'A': 1.0, 'x0': 0.5}, + parameter_units={'A': '1/meV', 'x0': 'meV'}, + x_unit='meV', + ) + + # THEN: relabelling A to 1/ueV makes the output meV/ueV, which converts to dimensionless + with warnings.catch_warnings(): + warnings.simplefilter('error') expr.set_unit('A', '1/ueV') + # EXPECT: evaluated values carry the 1000x conversion into y_unit + assert expr.evaluate(np.array([1.5]))[0] == pytest.approx(1000.0) + + def test_convertible_output_unit_is_rescaled_into_y_unit(self): + # WHEN: the jump-diffusion width in SI-flavoured parameter units, wanted in meV + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent( + 'hbar * D * x**2 / (1 + D * x**2 * tau)', + parameters={'D': 1e-9, 'tau': 1.0}, + parameter_units={'D': 'm^2/s', 'tau': 'ps'}, + x_unit='1/angstrom', + y_unit='meV', + ) + + # THEN + value = expr.evaluate(np.array([1.0]))[0] + + # EXPECT: hbar * D * Q^2 / (1 + D * Q^2 * tau) expressed in meV. With + # hbar = 6.582120e-13 meV*s, D = 1e-9 m^2/s = 1e11 angstrom^2/s and tau = 1e-12 s the + # denominator is 1 + 0.1 and the numerator 6.582120e-2 meV. + assert value == pytest.approx(6.582120e-2 / 1.1, rel=1e-5) + + def test_conversion_handles_non_si_dimensions_like_counts(self): + # WHEN: an intensity-scaled jump-diffusion width, wanted in counts*meV + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent( + 'counts * hbar * D * x**2 / (1 + D * x**2 * tau)', + parameters={'counts': 1.0, 'D': 4.6e-10, 'tau': 22.0}, + parameter_units={'counts': 'counts', 'D': 'm^2/s', 'tau': 'ps'}, + x_unit='1/angstrom', + y_unit='counts*meV', + ) + + # THEN + value = expr.evaluate(np.array([1.0]))[0] + + # EXPECT: counts is a non-SI dimension scipp carries in the unit powers; only the scale + # multiplier is converted. D*Q^2*tau = 4.6e-10 m^2/s * 1e20 /m^2 * 22e-12 s = 1.012. + expected = 6.582120e-13 * 4.6e-10 * 1e20 / (1.0 + 4.6e-10 * 1e20 * 22e-12) + assert value == pytest.approx(expected, rel=1e-5) + + ############# + # Physical constants + ############# -class TestExpressionComponentPhysicalConstants: def test_kb_constant_value_and_unit(self): # WHEN expr = ExpressionComponent( diff --git a/tests/unit/easydynamics/sample_model/components/test_gaussian.py b/tests/unit/easydynamics/sample_model/components/test_gaussian.py index 3ef9b01b7..266c6802c 100644 --- a/tests/unit/easydynamics/sample_model/components/test_gaussian.py +++ b/tests/unit/easydynamics/sample_model/components/test_gaussian.py @@ -131,6 +131,35 @@ def test_width_must_be_positive(self, gaussian: Gaussian): with pytest.raises(ValueError, match='width must be positive'): gaussian.width = -0.5 + def test_area_setter_out_of_bounds_raises(self, gaussian: Gaussian): + # WHEN the fixture's area was created non-negative, so it carries min=0 + original_area = gaussian.area.value + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + gaussian.area = -1.0 + assert gaussian.area.value == pytest.approx(original_area) + + def test_area_setter_allows_negative_when_unbounded(self): + # WHEN a Gaussian constructed with a negative area gets no lower bound + with pytest.warns(UserWarning, match='may not be physically meaningful'): + gaussian = Gaussian(area=-2.0) + + # THEN + gaussian.area = -1.0 + + # EXPECT + assert gaussian.area.value == pytest.approx(-1.0) + + def test_width_setter_below_minimum_raises(self, gaussian: Gaussian): + # WHEN the width parameter carries an absolute minimum (1e-10) + original_width = gaussian.width.value + + # THEN EXPECT a tiny positive width below the bound raises instead of being clamped + with pytest.raises(ValueError, match='violates the parameter bounds'): + gaussian.width = 1e-12 + assert gaussian.width.value == pytest.approx(original_width) + def test_evaluate(self, gaussian: Gaussian): # WHEN x = np.array([0.0, 0.5, 1.0]) diff --git a/tests/unit/easydynamics/sample_model/components/test_lorentzian.py b/tests/unit/easydynamics/sample_model/components/test_lorentzian.py index 97e02aad9..3a25d3073 100644 --- a/tests/unit/easydynamics/sample_model/components/test_lorentzian.py +++ b/tests/unit/easydynamics/sample_model/components/test_lorentzian.py @@ -124,6 +124,24 @@ def test_width_must_be_positive(self, lorentzian: Lorentzian): with pytest.raises(ValueError, match='width must be positive'): lorentzian.width = -0.5 + def test_area_setter_out_of_bounds_raises(self, lorentzian: Lorentzian): + # WHEN the fixture's area was created non-negative, so it carries min=0 + original_area = lorentzian.area.value + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + lorentzian.area = -1.0 + assert lorentzian.area.value == pytest.approx(original_area) + + def test_width_setter_below_minimum_raises(self, lorentzian: Lorentzian): + # WHEN the width parameter carries an absolute minimum (1e-10) + original_width = lorentzian.width.value + + # THEN EXPECT a tiny positive width below the bound raises instead of being clamped + with pytest.raises(ValueError, match='violates the parameter bounds'): + lorentzian.width = 1e-12 + assert lorentzian.width.value == pytest.approx(original_width) + def test_evaluate(self, lorentzian: Lorentzian): # WHEN x = np.array([0.0, 0.5, 1.0]) diff --git a/tests/unit/easydynamics/sample_model/components/test_mixins.py b/tests/unit/easydynamics/sample_model/components/test_mixins.py index d7a8487ca..d275e6666 100644 --- a/tests/unit/easydynamics/sample_model/components/test_mixins.py +++ b/tests/unit/easydynamics/sample_model/components/test_mixins.py @@ -13,7 +13,9 @@ class TestCreateParametersMixin: def dummy_model(self): return CreateParametersMixin() - # ------------- Area---------------------- + ############# + # Area + ############# @pytest.mark.parametrize('unit', ['meV', 'eV']) @pytest.mark.parametrize('area_input', [2, 2.0]) def test_create_area_parameter_from_numeric(self, dummy_model, area_input, unit): @@ -53,7 +55,40 @@ def test_negative_area_warns(self, dummy_model): assert area_param.min == -float('inf') # No min constraint for negative area - # ------------- Center---------------------- + ############# + # Bounded value assignment + ############# + def test_set_bounded_parameter_value_within_bounds(self, dummy_model): + # WHEN + param = Parameter(name='p', value=1.0, min=0.0, max=2.0) + + # THEN + dummy_model._set_bounded_parameter_value(param, 1.5, 'p') + + # EXPECT + assert param.value == pytest.approx(1.5) + + @pytest.mark.parametrize('out_of_bounds', [-1.0, 3.0], ids=['below_min', 'above_max']) + def test_set_bounded_parameter_value_out_of_bounds_raises(self, dummy_model, out_of_bounds): + # WHEN + param = Parameter(name='p', value=1.0, min=0.0, max=2.0) + + # THEN EXPECT the assignment raises instead of silently clamping, leaving the value + with pytest.raises(ValueError, match='violates the parameter bounds'): + dummy_model._set_bounded_parameter_value(param, out_of_bounds, 'p') + assert param.value == pytest.approx(1.0) + + def test_set_bounded_parameter_value_invalid_type_raises(self, dummy_model): + # WHEN + param = Parameter(name='p', value=1.0, min=0.0, max=2.0) + + # THEN EXPECT + with pytest.raises(TypeError, match='p must be a number'): + dummy_model._set_bounded_parameter_value(param, 'invalid', 'p') + + ############# + # Center + ############# @pytest.mark.parametrize('unit', ['meV', 'eV']) @pytest.mark.parametrize('center_input', [0, 0.0]) def test_create_center_parameter_from_numeric(self, dummy_model, center_input, unit): diff --git a/tests/unit/easydynamics/sample_model/components/test_model_component.py b/tests/unit/easydynamics/sample_model/components/test_model_component.py index 40e05e03a..f806afc75 100644 --- a/tests/unit/easydynamics/sample_model/components/test_model_component.py +++ b/tests/unit/easydynamics/sample_model/components/test_model_component.py @@ -231,7 +231,9 @@ def test_evaluate_with_compatible_unit_gives_correct_result(self): assert g_mev.width.value == pytest.approx(0.5) assert g_mev.area.value == pytest.approx(1.0) - # ───── Regression tests ───── + ############# + # Regression tests + ############# def test_convert_x_unit_rollback_on_failure(self, dummy: DummyComponent): # Conversion to 'm' (length) is incompatible with 'meV' (energy) → triggers rollback diff --git a/tests/unit/easydynamics/sample_model/components/test_polynomial.py b/tests/unit/easydynamics/sample_model/components/test_polynomial.py index a5ca6a599..692c32993 100644 --- a/tests/unit/easydynamics/sample_model/components/test_polynomial.py +++ b/tests/unit/easydynamics/sample_model/components/test_polynomial.py @@ -184,6 +184,34 @@ def test_convert_x_unit_raises_invalid_unit(self, polynomial: Polynomial): with pytest.raises(Exception, match='unit must be '): polynomial.convert_x_unit(123) + def test_convert_x_unit_rescales_bounded_coefficient_without_clamping(self): + # WHEN a linear coefficient with a lower bound that the converted value would cross + # (regression: the value was multiplied in place and easyscience silently clamped + # it to the bound, corrupting the coefficient irreversibly) + bounded = Parameter(name='c1', value=1.0, min=0.5) + polynomial = Polynomial(coefficients=[0.0, bounded], x_unit='meV') + + # THEN + polynomial.convert_x_unit('microeV') + + # EXPECT the value and the bound are rescaled together instead of clamping + assert bounded.value == pytest.approx(1e-3) + assert bounded.min == pytest.approx(0.5e-3) + # and the evaluated polynomial is physically unchanged: 1000 microeV = 1 meV + assert polynomial.evaluate(np.array([1000.0]))[0] == pytest.approx(1.0) + + def test_convert_y_unit_rescales_bounded_coefficient_without_clamping(self): + # WHEN a coefficient with an upper bound that the converted value would cross + bounded = Parameter(name='c0', value=1.0, max=2.0) + polynomial = Polynomial(coefficients=[bounded], x_unit='meV', y_unit='1/meV') + + # THEN + polynomial.convert_y_unit('1/eV') + + # EXPECT the value and the bound are rescaled together instead of clamping + assert bounded.value == pytest.approx(1e3) + assert bounded.max == pytest.approx(2e3) + def test_copy(self, polynomial: Polynomial): # WHEN THEN polynomial_copy = copy(polynomial) @@ -299,7 +327,9 @@ def test_convert_y_unit_rollback_on_failure(self): assert np.isclose(p.coefficients[0].value, 1.0) assert np.isclose(p.coefficients[1].value, 2.0) - # --- Serialization --- + ############# + # Serialization + ############# def test_to_dict(self, polynomial: Polynomial): # WHEN @@ -359,7 +389,9 @@ def test_from_dict_invalid_dict_raises(self): with pytest.raises(ValueError, match='must be a dictionary representing'): Polynomial.from_dict({'not': 'valid'}) - # --- Sparse dict initialization --- + ############# + # Sparse dict initialization + ############# def test_sparse_dict_single_term(self): # WHEN @@ -421,7 +453,9 @@ def test_sparse_dict_preserves_x_unit(self): assert p.x_unit == 'ueV' assert p.y_unit == 'counts' - # --- add_coefficient --- + ############# + # add_coefficient + ############# def test_add_coefficient_increases_degree(self, polynomial: Polynomial): # WHEN @@ -470,7 +504,9 @@ def test_add_coefficient_appears_in_all_variables(self, polynomial: Polynomial): # THEN EXPECT: the new coefficient is fittable assert len(polynomial.get_all_variables()) == 4 - # --- remove_coefficient --- + ############# + # remove_coefficient + ############# def test_remove_coefficient_decreases_degree(self, polynomial: Polynomial): # WHEN @@ -533,10 +569,9 @@ def test_add_coefficient_then_convert_x_unit(self, polynomial: Polynomial): assert polynomial.x_unit == 'ueV' np.testing.assert_allclose(after, before, rtol=1e-8) - -def test_suppress_warnings_setter_raises_for_non_bool(): - # GIVEN a Polynomial - p = Polynomial(coefficients=[1.0, 2.0], x_unit='meV') - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'Suppress_warnings must be True or False'): - p.suppress_warnings = 'yes' + def test_suppress_warnings_setter_raises_for_non_bool(self): + # WHEN a Polynomial + p = Polynomial(coefficients=[1.0, 2.0], x_unit='meV') + # THEN EXPECT + with pytest.raises(TypeError, match=r'Suppress_warnings must be True or False'): + p.suppress_warnings = 'yes' diff --git a/tests/unit/easydynamics/sample_model/components/test_voigt.py b/tests/unit/easydynamics/sample_model/components/test_voigt.py index 42eb5099c..8c66b886e 100644 --- a/tests/unit/easydynamics/sample_model/components/test_voigt.py +++ b/tests/unit/easydynamics/sample_model/components/test_voigt.py @@ -201,6 +201,24 @@ def test_lorentzian_width_must_be_positive(self, voigt: Voigt): ): voigt.lorentzian_width = -0.7 + def test_area_setter_out_of_bounds_raises(self, voigt: Voigt): + # WHEN the fixture's area was created non-negative, so it carries min=0 + original_area = voigt.area.value + + # THEN EXPECT a negative assignment raises instead of being silently clamped to 0 + with pytest.raises(ValueError, match='violates the parameter bounds'): + voigt.area = -1.0 + assert voigt.area.value == pytest.approx(original_area) + + def test_width_setters_below_minimum_raise(self, voigt: Voigt): + # WHEN the width parameters carry an absolute minimum (1e-10) + + # THEN EXPECT tiny positive widths below the bound raise instead of being clamped + with pytest.raises(ValueError, match='violates the parameter bounds'): + voigt.gaussian_width = 1e-12 + with pytest.raises(ValueError, match='violates the parameter bounds'): + voigt.lorentzian_width = 1e-12 + def test_center_is_fixed_if_set_to_None(self, voigt: Voigt): # WHEN assert voigt.center.fixed is False diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_brownian_translational_diffusion.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_brownian_translational_diffusion.py index 4b2137540..c8958afbb 100644 --- a/tests/unit/easydynamics/sample_model/diffusion_model/test_brownian_translational_diffusion.py +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_brownian_translational_diffusion.py @@ -4,7 +4,6 @@ import numpy as np import pytest import scipp as sc -from easyscience.variable import DescriptorNumber from scipp import UnitError from scipp.constants import hbar as scipp_hbar @@ -12,10 +11,6 @@ BrownianTranslationalDiffusion, ) -hbar_1 = DescriptorNumber('hbar', 1.0) -hbar = DescriptorNumber.from_scipp('hbar', scipp_hbar) -angstrom = DescriptorNumber('angstrom', 1e-10, unit='m') - class TestBrownianTranslationalDiffusion: @pytest.fixture @@ -258,6 +253,17 @@ def test_create_component_collections(self, brownian_diffusion_model, Q): # area.unit = area_unit = x_unit * y_unit assert component.area.unit == 'meV' + def test_create_component_collections_installs_collections(self): + # WHEN + model = BrownianTranslationalDiffusion(Q=np.array([1.0, 2.0])) + + # THEN + collections = model.create_component_collections() + + # EXPECT the returned collections are the installed (live) ones, so callers that + # follow the docstring get the same objects the model itself uses + assert collections is model.get_component_collections() + def test_write_width_dependency_expression(self, brownian_diffusion_model): # WHEN THEN expression = brownian_diffusion_model._write_width_dependency_expression(0.5) diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_delta_lorentz.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_delta_lorentz.py index 17074eb95..59f2d03f9 100644 --- a/tests/unit/easydynamics/sample_model/diffusion_model/test_delta_lorentz.py +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_delta_lorentz.py @@ -265,9 +265,9 @@ def test_input_type_validation_raises(self, kwargs, expected_exception, expected with pytest.raises(expected_exception, match=expected_message): DeltaLorentz(**kwargs) - # ------------------------------------------------------------------ + ############# # Properties - # ------------------------------------------------------------------ + ############# @pytest.mark.parametrize( ('attribute', 'value', 'expected'), [ @@ -418,9 +418,9 @@ def test_setters_invalid( with pytest.raises(exception, match=message): setattr(delta_lorentz_model, attribute, value) - # ------------------------------------------------------------------ + ############# # Other methods - # ------------------------------------------------------------------ + ############# def test_calculate_width_without_Q(self, delta_lorentz_model): # WHEN THEN @@ -565,6 +565,47 @@ def test_create_component_collections_with_no_Q_variation( ) assert 'A_0' in collection[1].area.dependency_expression + def test_create_component_collections_installs_and_stays_in_sync( + self, delta_lorentz_model_with_Q + ): + # WHEN + model = delta_lorentz_model_with_Q + + # THEN + collections = model.create_component_collections() + + # EXPECT the returned collections are the installed (live) ones (regression: they + # were returned without being installed, while the per-Q parameter lists were + # replaced, desynchronizing calculate_width from the installed components) + assert collections is model.get_component_collections() + + # THEN setting a per-Q width parameter + model._lorentzian_width_list[0].value = 0.5 + + # EXPECT the change is visible in the installed component and in calculate_width + assert collections[0][0].width.value == pytest.approx(0.5) + assert model.calculate_width()[0] == pytest.approx(0.5) + + # EXPECT the same holds for the per-Q amplitude parameters + model._A_0_list[0].value = 0.25 + assert model.calculate_EISF()[0] == pytest.approx(0.25) + assert collections[0][1].area.value == pytest.approx(0.25) + + def test_per_Q_parameter_and_collection_names(self, delta_lorentz_model_with_Q): + # WHEN + model = delta_lorentz_model_with_Q + + # THEN + collections = model.get_component_collections() + + # EXPECT the per-Q amplitudes carry the model name (like the widths carry the + # Lorentzian name), and the per-Q collections get a name, not just a display name + for a0, a1 in zip(model._A_0_list, model._A_1_list, strict=True): + assert a0.name == 'DeltaLorentz A_0' + assert a1.name == 'DeltaLorentz A_1' + assert collections[0].name == 'DeltaLorentz_Q0.50' + assert collections[0].display_name == 'DeltaLorentz_Q0.50' + @pytest.mark.parametrize( ('Q_index', 'expected_exception', 'expected_message'), [ @@ -923,7 +964,9 @@ def test_repr(self, delta_lorentz_model): # Regression: a stray ')' used to mangle this into 'x_unit=meV), y_unit=...' assert 'x_unit=meV, y_unit=dimensionless' in repr_str - # ───── Regression tests ───── + ############# + # Regression tests + ############# def test_calculate_width_with_Q_subset(self, delta_lorentz_model_with_Q): # WHEN: Q-varying widths with distinguishable per-Q values @@ -972,21 +1015,23 @@ def test_calculate_width_raises_after_clear_Q_when_allow_Q_variation( with pytest.raises(ValueError, match='Q must be provided'): delta_lorentz_model_with_Q.calculate_width() + ############# + # Fit targets and Q validation + ############# -def test_get_fit_targets_includes_delta_area(): - # GIVEN a DeltaLorentz model - model = DeltaLorentz(delta_name='Delta function', lorentzian_name='Lorentzian') - # WHEN - targets = model.get_fit_targets() - # EXPECT base area/width plus the delta_area prediction - assert [t.name for t in targets] == ['area', 'width', 'delta_area'] - delta_area = next(t for t in targets if t.name == 'delta_area') - assert delta_area.dataset_key == 'Delta function area' - - -def test_calculate_width_raises_when_Q_variation_enabled_but_Q_unset(): - # GIVEN Q-variation enabled for the width but Q never set on the model (empty per-Q list) - model = DeltaLorentz(lorentzian_width=0.1, allow_Q_variation={'lorentzian_width': True}) - # WHEN a Q is requested THEN EXPECT the empty per-Q width list to be reported - with pytest.raises(ValueError, match=r'Lorentzian width Q-variation list is empty'): - model.calculate_width(np.array([1.0])) + def test_get_fit_targets_includes_delta_area(self): + # WHEN a DeltaLorentz model + model = DeltaLorentz(delta_name='Delta function', lorentzian_name='Lorentzian') + # THEN + targets = model.get_fit_targets() + # EXPECT base area/width plus the delta_area prediction + assert [t.name for t in targets] == ['area', 'width', 'delta_area'] + delta_area = next(t for t in targets if t.name == 'delta_area') + assert delta_area.dataset_key == 'Delta function area' + + def test_calculate_width_raises_when_Q_variation_enabled_but_Q_unset(self): + # WHEN Q-variation enabled for the width but Q never set on the model (empty per-Q list) + model = DeltaLorentz(lorentzian_width=0.1, allow_Q_variation={'lorentzian_width': True}) + # THEN EXPECT the empty per-Q width list to be reported when a Q is requested + with pytest.raises(ValueError, match=r'Lorentzian width Q-variation list is empty'): + model.calculate_width(np.array([1.0])) diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model_base.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model_base.py index d323b0ba6..ace30385e 100644 --- a/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model_base.py +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model_base.py @@ -355,21 +355,19 @@ def test_ensure_Q_uses_argument(self, diffusion_model): # EXPECT np.testing.assert_allclose(Q, [1.0, 2.0]) - -def test_get_fit_targets_declares_area_and_width(): - # GIVEN a diffusion model - model = DiffusionModelBase(lorentzian_name='Lorentzian') - # WHEN - targets = model.get_fit_targets() - # EXPECT area and width predictions with keys derived from the Lorentzian name - assert [t.name for t in targets] == ['area', 'width'] - assert targets[0].dataset_key == 'Lorentzian area' - assert targets[1].dataset_key == 'Lorentzian width' - - -def test_match_Q_indices_raises_when_Q_not_set(): - # GIVEN a diffusion model with no Q set - model = DiffusionModelBase() - # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'Q must be set in the model'): - model._match_Q_indices(np.array([1.0])) + def test_get_fit_targets_declares_area_and_width(self): + # WHEN a diffusion model + model = DiffusionModelBase(lorentzian_name='Lorentzian') + # THEN + targets = model.get_fit_targets() + # EXPECT area and width predictions with keys derived from the Lorentzian name + assert [t.name for t in targets] == ['area', 'width'] + assert targets[0].dataset_key == 'Lorentzian area' + assert targets[1].dataset_key == 'Lorentzian width' + + def test_match_Q_indices_raises_when_Q_not_set(self): + # WHEN a diffusion model with no Q set + model = DiffusionModelBase() + # THEN EXPECT + with pytest.raises(ValueError, match=r'Q must be set in the model'): + model._match_Q_indices(np.array([1.0])) diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_jump_translational_diffusion.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_jump_translational_diffusion.py index 7fcf4dacd..4aa53200b 100644 --- a/tests/unit/easydynamics/sample_model/diffusion_model/test_jump_translational_diffusion.py +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_jump_translational_diffusion.py @@ -4,7 +4,6 @@ import numpy as np import pytest import scipp as sc -from easyscience.variable import DescriptorNumber from scipp import UnitError from scipp.constants import hbar as scipp_hbar @@ -12,10 +11,6 @@ JumpTranslationalDiffusion, ) -hbar_1 = DescriptorNumber('hbar', 1.0) -hbar = DescriptorNumber.from_scipp('hbar', scipp_hbar) -angstrom = DescriptorNumber('angstrom', 1e-10, unit='m') - class TestJumpTranslationalDiffusion: @pytest.fixture @@ -248,6 +243,17 @@ def test_create_component_collections(self, jump_diffusion_model, Q): # area.unit = area_unit = x_unit * y_unit assert component.area.unit == 'meV' + def test_create_component_collections_installs_collections(self): + # WHEN + model = JumpTranslationalDiffusion(Q=np.array([1.0, 2.0])) + + # THEN + collections = model.create_component_collections() + + # EXPECT the returned collections are the installed (live) ones, so callers that + # follow the docstring get the same objects the model itself uses + assert collections is model.get_component_collections() + def test_write_width_dependency_expression(self, jump_diffusion_model): # WHEN THEN expression = jump_diffusion_model._write_width_dependency_expression(0.5) diff --git a/tests/unit/easydynamics/sample_model/test_component_collection.py b/tests/unit/easydynamics/sample_model/test_component_collection.py index 233117031..be3407f14 100644 --- a/tests/unit/easydynamics/sample_model/test_component_collection.py +++ b/tests/unit/easydynamics/sample_model/test_component_collection.py @@ -9,6 +9,7 @@ from easyscience.variable import Parameter from scipy.integrate import simpson +from easydynamics.exceptions import AmbiguousNameError from easydynamics.sample_model import ComponentCollection from easydynamics.sample_model import ExpressionComponent from easydynamics.sample_model import Gaussian @@ -113,7 +114,9 @@ def test_init_with_invalid_unit_raises(self): with pytest.raises(TypeError, match='unit must be'): ComponentCollection(x_unit=123) - # ───── Component Management ───── + ############# + # Component Management + ############# def test_append_component(self, component_collection): # WHEN @@ -308,7 +311,9 @@ def test_evaluate_component_invalid_name_type_raises(self, component_collection) ): component_collection.evaluate_component(x, 123) - # ───── Utilities ───── + ############# + # Utilities + ############# def test_normalize_area(self, component_collection): # WHEN THEN @@ -676,7 +681,95 @@ def test_evaluate_scipp_output_with_y_unit(self): assert isinstance(result, sc.Variable) assert result.unit == sc.Unit('1/meV') - # ───── Regression tests ───── + ############# + # Versioning + ############# + + def test_version_starts_at_zero_and_bumps_on_mutation(self): + # WHEN a freshly constructed collection with initial components + collection = ComponentCollection(components=[Gaussian(name='G1'), Lorentzian(name='L1')]) + + # EXPECT it starts at version 0 + assert collection.version == 0 + + # THEN structural mutations bump the version + collection.append_component(Gaussian(name='G2')) + assert collection.version == 1 + collection.pop('G2') + assert collection.version == 2 + + ############# + # Slicing + ############# + + def test_getitem_slice_returns_working_collection(self, component_collection): + "Regression: slicing used to crash because the base slice path called the wrong ctor" + # WHEN THEN + sliced = component_collection[:1] + + # EXPECT a working collection of the same class, carrying the units, sharing the + # component objects + assert type(sliced) is ComponentCollection + assert len(sliced) == 1 + assert sliced[0] is component_collection[0] + assert sliced.x_unit == component_collection.x_unit + assert sliced.y_unit == component_collection.y_unit + + # EXPECT the slice is usable + x = np.linspace(-5, 5, 11) + np.testing.assert_allclose(sliced.evaluate(x), component_collection[0].evaluate(x)) + + ############# + # Regression tests + ############# + + def test_normalize_area_negative_area_raises(self, component_collection): + "Regression: negative areas used to be silently clamped by normalization" + # WHEN + component_collection[0].area.min = -10.0 + component_collection[0].area = -2.0 + + # THEN EXPECT + with pytest.raises(ValueError, match=r'Negative area'): + component_collection.normalize_area() + + def test_evaluate_empty_invalid_output_raises(self): + "Regression: the empty-collection path used to skip output validation" + # WHEN + collection = ComponentCollection(display_name='EmptyModel') + + # THEN EXPECT + with pytest.raises(ValueError, match=r"output must be 'numpy' or 'scipp'"): + collection.evaluate(np.linspace(-1, 1, 5), output='invalid') + + def test_evaluate_empty_scalar_shape_matches_non_empty_path(self): + "Regression: empty and non-empty paths must agree on the output shape for scalar x" + # WHEN an empty and a non-empty collection evaluated at a scalar + empty = ComponentCollection(display_name='EmptyModel') + non_empty = ComponentCollection(components=Gaussian(name='G')) + + # THEN + empty_result = empty.evaluate(0.5) + non_empty_result = non_empty.evaluate(0.5) + + # EXPECT both return 1D arrays of the same shape + assert empty_result.shape == non_empty_result.shape == (1,) + assert np.all(empty_result == pytest.approx(0.0)) + + def test_evaluate_component_ambiguous_name_raises(self): + "Regression: duplicate names used to silently evaluate the first match" + # WHEN a collection with two components sharing a name + with pytest.warns(UserWarning, match='Duplicate component names'): + collection = ComponentCollection( + components=[ + Gaussian(name='SameName', area=1.0), + Gaussian(name='SameName', area=2.0), + ] + ) + + # THEN EXPECT + with pytest.raises(AmbiguousNameError, match=r"Ambiguous name 'SameName'"): + collection.evaluate_component(np.linspace(-1, 1, 5), 'SameName') def test_evaluate_scipp_output_multi_component_does_not_raise(self, component_collection): # WHEN: collection with two components (Gaussian + Lorentzian) diff --git a/tests/unit/easydynamics/sample_model/test_instrument_model.py b/tests/unit/easydynamics/sample_model/test_instrument_model.py index a802b1961..6df773e10 100644 --- a/tests/unit/easydynamics/sample_model/test_instrument_model.py +++ b/tests/unit/easydynamics/sample_model/test_instrument_model.py @@ -551,36 +551,39 @@ def test_on_energy_offset_change(self, instrument_model): assert offset.value == new_offset def test_on_resolution_model_change(self, instrument_model, resolution_model): - # WHEN + # WHEN a resolution model that does not know the instrument's Q yet new_resolution_model = resolution_model + assert new_resolution_model.Q is None # THEN - instrument_model._resolution_model = new_resolution_model - instrument_model._on_resolution_model_change() + instrument_model.resolution_model = new_resolution_model - # EXPECT - assert instrument_model._resolution_model is new_resolution_model + # EXPECT the change callback propagated the instrument's Q to the new model + assert instrument_model.resolution_model is new_resolution_model + np.testing.assert_array_equal(new_resolution_model.Q.values, np.array([1.0, 2.0, 3.0])) def test_on_background_model_change(self, instrument_model, background_model): - # WHEN + # WHEN a background model that does not know the instrument's Q yet new_background_model = background_model + assert new_background_model.Q is None # THEN - instrument_model._background_model = new_background_model - instrument_model._on_background_model_change() + instrument_model.background_model = new_background_model - # EXPECT - assert instrument_model._background_model is new_background_model + # EXPECT the change callback propagated the instrument's Q to the new model + assert instrument_model.background_model is new_background_model + np.testing.assert_array_equal(new_background_model.Q.values, np.array([1.0, 2.0, 3.0])) def test_repr_contains_expected_fields(self, instrument_model): # WHEN THEN repr_str = repr(instrument_model) - # EXPECT + # EXPECT values pinned from the fixture's known construction inputs, so a wrong + # attribute value cannot satisfy its own interpolation assert repr_str.startswith('InstrumentModel(') - assert f'unique_name={instrument_model.unique_name!r}' in repr_str - assert f'x_unit={instrument_model.x_unit}' in repr_str - assert 'Q_len=3' in repr_str - assert f'resolution_model={instrument_model._resolution_model!r}' in repr_str - assert f'background_model={instrument_model._background_model!r}' in repr_str assert repr_str.endswith(')') + assert "unique_name='" in repr_str + assert 'x_unit=meV' in repr_str + assert 'Q_len=3' in repr_str + assert 'resolution_model=ResolutionModel(' in repr_str + assert 'background_model=BackgroundModel(' in repr_str diff --git a/tests/unit/easydynamics/sample_model/test_model_base.py b/tests/unit/easydynamics/sample_model/test_model_base.py index 670895edb..5bd222185 100644 --- a/tests/unit/easydynamics/sample_model/test_model_base.py +++ b/tests/unit/easydynamics/sample_model/test_model_base.py @@ -543,3 +543,76 @@ def test_convert_y_unit_invalid_raises(self, model_base): # WHEN THEN EXPECT with pytest.raises(TypeError): model_base.convert_y_unit(123) + + ############# + # State versioning + ############# + + def test_evaluate_without_Q_names_the_cause(self): + "Regression: the error used to claim 'no components' when Q was the missing piece" + # WHEN a model with components but no Q + model = ModelBase(display_name='M', components=Gaussian(name='G')) + + # THEN EXPECT + with pytest.raises(ValueError, match='Q is not set'): + model.evaluate(np.array([0.0, 1.0])) + + def test_state_version_reading_does_not_mutate(self, model_base): + # WHEN + version = model_base.state_version + + # THEN EXPECT repeated reads return the same value and rebuild nothing + assert model_base.state_version == version + assert model_base.component_collections_is_dirty is True + assert model_base._component_collections == [] + + def test_state_version_changes_on_component_and_Q_changes(self, model_base): + # WHEN + version = model_base.state_version + + # THEN appending a component through the model + model_base.append_component(Gaussian(name='SVGaussian')) + # EXPECT + assert model_base.state_version > version + version = model_base.state_version + + # THEN removing a component through the model + model_base.remove_component('SVGaussian') + # EXPECT + assert model_base.state_version > version + version = model_base.state_version + + # THEN clearing Q + model_base.clear_Q(confirm=True) + # EXPECT + assert model_base.state_version > version + + def test_state_version_changes_on_in_place_template_mutation(self, model_base): + # WHEN collections are current, so the dirty flag alone would report clean + _ = model_base.get_component_collection(0) + assert model_base.component_collections_is_dirty is False + version = model_base.state_version + + # THEN mutating the live template collection in place, bypassing the model's methods + model_base.components.append_component(Gaussian(name='LiveGaussian')) + + # EXPECT the mutation is visible without any callback + assert model_base.state_version > version + assert model_base.component_collections_is_dirty is True + + def test_evaluate_includes_component_appended_to_live_collection(self, model_base): + "Regression: components appended via the live template collection were invisible" + # WHEN a model whose collections were already built and evaluated + x = np.linspace(-5, 5, 101) + result_before = model_base.evaluate(x) + + # THEN appending directly to the live template collection and evaluating again + model_base.components.append_component( + Gaussian(name='LiveGaussian', area=10.0, center=0.0, width=1.0) + ) + result_after = model_base.evaluate(x) + + # EXPECT the new component contributes to the output at every Q + extra = Gaussian(name='Reference', area=10.0, center=0.0, width=1.0).evaluate(x) + for before, after in zip(result_before, result_after, strict=True): + np.testing.assert_allclose(after, before + extra, rtol=1e-10) diff --git a/tests/unit/easydynamics/sample_model/test_resolution_model.py b/tests/unit/easydynamics/sample_model/test_resolution_model.py index 102363f60..5a5d2773d 100644 --- a/tests/unit/easydynamics/sample_model/test_resolution_model.py +++ b/tests/unit/easydynamics/sample_model/test_resolution_model.py @@ -305,18 +305,61 @@ def test_from_sample_model_invalid_arguments( **valid_kwargs, ) - def test_from_sample_model_invalid_components(self, sample_model): - # WHEN - invalid_component = DeltaFunction(name='InvalidDelta') - sample_model.append_component(invalid_component) + def test_from_sample_model_strips_delta_functions(self, sample_model): + # WHEN a sample model with the standard QENS elastic delta line + sample_model.append_component(DeltaFunction(name='Elastic')) - # THEN EXPECT - with pytest.raises( - TypeError, - match='cannot be a DeltaFunction', - ): + # THEN + with pytest.warns(UserWarning, match='Stripped'): + resolution_model = ResolutionModel.from_sample_model(sample_model) + + # EXPECT no DeltaFunction in the template or the per-Q collections, and the + # remaining components still normalized to unit area + assert not any(isinstance(c, DeltaFunction) for c in resolution_model.components) + for Q_index in range(len(resolution_model.Q)): + collection = resolution_model.get_component_collection(Q_index) + assert not any(isinstance(c, DeltaFunction) for c in collection) + assert sum(c.area.value for c in collection) == pytest.approx(1.0) + + def test_from_sample_model_background_component_raises(self, sample_model): + # WHEN a sample model carrying a background component + sample_model.append_component(Polynomial(name='Background')) + + # THEN EXPECT backgrounds are rejected, not silently installed as resolution + with pytest.raises(TypeError, match='cannot be a Polynomial'): + ResolutionModel.from_sample_model(sample_model) + + def test_from_sample_model_delta_only_raises(self): + # WHEN a sample model whose only component is the elastic delta + sample_model = SampleModel( + components=DeltaFunction(name='Elastic'), + Q=np.array([1.0]), + ) + + # THEN EXPECT stripping the delta would leave no resolution shape + with pytest.raises(ValueError, match='contains only'): ResolutionModel.from_sample_model(sample_model) + def test_from_sample_model_locks_calibrated_collections(self, sample_model): + # WHEN + resolution_model = ResolutionModel.from_sample_model(sample_model) + calibrated = resolution_model.get_component_collection(0) + + # THEN EXPECT mutations that would rebuild the collections from the unfitted + # template fail loudly instead of silently discarding the calibration + with pytest.raises(RuntimeError, match='calibrated'): + resolution_model.append_component(Gaussian(name='Extra')) + with pytest.raises(RuntimeError, match='calibrated'): + resolution_model.remove_component('TestGaussian1Name') + with pytest.raises(RuntimeError, match='calibrated'): + resolution_model.clear_components() + with pytest.raises(RuntimeError, match='calibrated'): + resolution_model.clear_Q(confirm=True) + + # EXPECT the calibrated collections survive untouched + assert resolution_model.get_component_collection(0) is calibrated + assert resolution_model._component_collections_is_dirty is False + def test_y_unit_setter_raises(self, resolution_model): # WHEN / THEN / EXPECT with pytest.raises(AttributeError): diff --git a/tests/unit/easydynamics/sample_model/test_sample_model.py b/tests/unit/easydynamics/sample_model/test_sample_model.py index 36ac3fb72..ab986e160 100644 --- a/tests/unit/easydynamics/sample_model/test_sample_model.py +++ b/tests/unit/easydynamics/sample_model/test_sample_model.py @@ -517,6 +517,82 @@ def test_evaluate_doesnt_call_dbf_when_disabled( np.testing.assert_allclose(result[0], np.array([1.0, 2.0, 3.0])) np.testing.assert_allclose(result[1], np.array([4.0, 5.0, 6.0])) + def test_evaluate_scipp_output_with_detailed_balance(self, sample_model): + # WHEN the fixture has temperature set, so detailed balance is applied + x = np.linspace(-2.0, 2.0, 21) + + # THEN (regression: multiplying an sc.Variable with the numpy DBF used to raise) + balanced = sample_model.evaluate(x, output='scipp') + + # EXPECT scipp output matches numpy output, with the model's y_unit kept + reference = sample_model.evaluate(x, output='numpy') + assert len(balanced) == 3 + for scipp_values, numpy_values in zip(balanced, reference, strict=True): + assert isinstance(scipp_values, sc.Variable) + assert scipp_values.unit == sc.Unit('dimensionless') + np.testing.assert_allclose(scipp_values.values, numpy_values) + + # THEN disabling detailed balance + sample_model.use_detailed_balance = False + unbalanced = sample_model.evaluate(x, output='scipp') + + # EXPECT the detailed balance factor really was applied above + assert not np.allclose(balanced[0].values, unbalanced[0].values) + + def test_evaluate_dataarray_input_with_and_without_detailed_balance(self, sample_model): + # WHEN + x = np.linspace(-2.0, 2.0, 21) + data_array = sc.DataArray( + data=sc.array(dims=['energy'], values=np.zeros_like(x)), + coords={'energy': sc.array(dims=['energy'], values=x, unit='meV')}, + ) + + # THEN (regression: detailed balance used to reject DataArray x, which the + # component pipeline explicitly supports) + with_temperature = sample_model.evaluate(data_array) + reference_with = sample_model.evaluate(x) + + sample_model.temperature = None + without_temperature = sample_model.evaluate(data_array) + reference_without = sample_model.evaluate(x) + + # EXPECT DataArray input matches plain numpy input in both modes + for result, reference in zip(with_temperature, reference_with, strict=True): + np.testing.assert_allclose(result, reference) + for result, reference in zip(without_temperature, reference_without, strict=True): + np.testing.assert_allclose(result, reference) + + def test_init_invalid_temperature_does_not_mutate_diffusion_models(self): + # WHEN a diffusion model without Q and an invalid temperature + diffusion_model = BrownianTranslationalDiffusion() + + # THEN EXPECT construction fails on the temperature validation + with pytest.raises(TypeError, match='temperature must be a number or None'): + SampleModel( + diffusion_models=diffusion_model, + Q=np.array([1.0, 2.0]), + temperature='cold', + ) + + # EXPECT the failed construction did not mutate the passed diffusion model + assert diffusion_model.Q is None + assert diffusion_model.get_component_collections() == [] + + def test_temperature_unit_is_normalized_to_str(self, sample_model): + # WHEN constructed with a scipp Unit instead of a string + model = SampleModel(temperature=10.0, temperature_unit=sc.Unit('K')) + + # EXPECT the stored unit is normalized to a string + assert isinstance(model.temperature_unit, str) + assert model.temperature_unit == 'K' + + # THEN converting with a scipp Unit + sample_model.convert_temperature_unit(sc.Unit('mK')) + + # EXPECT the stored unit is normalized to a string as well + assert isinstance(sample_model.temperature_unit, str) + assert sample_model.temperature_unit == 'mK' + def test_generate_component_collections(self, sample_model): # WHEN THEN sample_model._generate_component_collections() @@ -615,33 +691,31 @@ def test_convert_y_unit(self): assert model.components[0].y_unit == '1/eV' assert g.area.value == pytest.approx(1e3) - -def test_remove_diffusion_model_raises_with_duplicate_names(): - # GIVEN a SampleModel with two DiffusionModels sharing a name - Q = np.linspace(0.5, 2.0, 3) - model = SampleModel( - Q=Q, - diffusion_models=[ - BrownianTranslationalDiffusion(name='dup'), - BrownianTranslationalDiffusion(name='dup'), - ], - ) - # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'Multiple DiffusionModels share the name'): - model.remove_diffusion_model('dup') - - -def test_convert_x_unit_rolls_back_when_diffusion_model_conversion_fails(): - # GIVEN a SampleModel whose diffusion model raises during x-unit conversion - Q = np.linspace(0.5, 2.0, 3) - brownian = BrownianTranslationalDiffusion() - model = SampleModel(Q=Q, diffusion_models=brownian) - original_unit = model.x_unit - # WHEN the conversion fails partway through - with ( - patch.object(brownian, 'convert_x_unit', side_effect=RuntimeError('boom')), - pytest.raises(RuntimeError, match='boom'), - ): - model.convert_x_unit('ueV') - # EXPECT the model's own x_unit to be rolled back to the original - assert model.x_unit == original_unit + def test_remove_diffusion_model_raises_with_duplicate_names(self): + # WHEN a SampleModel with two DiffusionModels sharing a name + Q = np.linspace(0.5, 2.0, 3) + model = SampleModel( + Q=Q, + diffusion_models=[ + BrownianTranslationalDiffusion(name='dup'), + BrownianTranslationalDiffusion(name='dup'), + ], + ) + # THEN EXPECT + with pytest.raises(ValueError, match=r'Multiple DiffusionModels share the name'): + model.remove_diffusion_model('dup') + + def test_convert_x_unit_rolls_back_when_diffusion_model_conversion_fails(self): + # WHEN a SampleModel whose diffusion model raises during x-unit conversion + Q = np.linspace(0.5, 2.0, 3) + brownian = BrownianTranslationalDiffusion() + model = SampleModel(Q=Q, diffusion_models=brownian) + original_unit = model.x_unit + # THEN EXPECT the conversion fails partway through + with ( + patch.object(brownian, 'convert_x_unit', side_effect=RuntimeError('boom')), + pytest.raises(RuntimeError, match='boom'), + ): + model.convert_x_unit('ueV') + # EXPECT the model's own x_unit to be rolled back to the original + assert model.x_unit == original_unit diff --git a/tests/unit/easydynamics/settings/test_convolution_settings.py b/tests/unit/easydynamics/settings/test_convolution_settings.py index a5c0ceb38..d5d210202 100644 --- a/tests/unit/easydynamics/settings/test_convolution_settings.py +++ b/tests/unit/easydynamics/settings/test_convolution_settings.py @@ -168,16 +168,25 @@ def test_extension_factor_setter_valid(self, default_convolution_settings, value assert default_convolution_settings.extension_factor == pytest.approx(float(value)) assert default_convolution_settings._plan_version == version_before + 1 + def test_extension_factor_setter_none(self, default_convolution_settings): + # WHEN + version_before = default_convolution_settings._plan_version + + # THEN None is accepted, matching the None-capable constructor + default_convolution_settings.extension_factor = None + + # EXPECT: value stored and the plan invalidated for all convolvers + assert default_convolution_settings.extension_factor is None + assert default_convolution_settings._plan_version == version_before + 1 + @pytest.mark.parametrize( 'value, expected_exception, match', [ ('0.2', TypeError, 'must be a number'), - (None, TypeError, 'must be a number'), (-0.1, ValueError, 'must be non-negative'), ], ids=[ 'not_numeric', - 'none', 'negative', ], ) diff --git a/tests/unit/easydynamics/settings/test_detailed_balance_settings.py b/tests/unit/easydynamics/settings/test_detailed_balance_settings.py index dba4d9d68..fa7821d35 100644 --- a/tests/unit/easydynamics/settings/test_detailed_balance_settings.py +++ b/tests/unit/easydynamics/settings/test_detailed_balance_settings.py @@ -109,6 +109,36 @@ def test_setters_invalid( with pytest.raises(expected_exception, match=match): default_detailed_balance_settings.normalize_detailed_balance = value + ############# + # Plan invalidation + ############# + + def test_setters_bump_plan_version(self, default_detailed_balance_settings): + "Regression: flag toggles used to be invisible to convolvers holding these settings" + # WHEN + settings = default_detailed_balance_settings + version_before = settings._plan_version + + # THEN toggling each flag + settings.use_detailed_balance = False + settings.normalize_detailed_balance = False + + # EXPECT one bump per changed flag + assert settings._plan_version == version_before + 2 + + def test_invalidate_plan_bumps_version(self, default_detailed_balance_settings): + # WHEN + settings = default_detailed_balance_settings + version_before = settings._plan_version + + # THEN + settings._invalidate_plan() + + # EXPECT + assert settings._plan_version == version_before + 1 + assert settings._plan_valid_for(version_before) is False + assert settings._plan_valid_for(settings._plan_version) is True + def test_repr_default(self, default_detailed_balance_settings): # WHEN repr_str = repr(default_detailed_balance_settings) diff --git a/tests/unit/easydynamics/test_exceptions.py b/tests/unit/easydynamics/test_exceptions.py index bc7731a5e..8547024c6 100644 --- a/tests/unit/easydynamics/test_exceptions.py +++ b/tests/unit/easydynamics/test_exceptions.py @@ -3,6 +3,7 @@ from easydynamics.exceptions import AmbiguousNameError +from easydynamics.sample_model import Gaussian class TestAmbiguousNameError: @@ -21,6 +22,22 @@ def test_initialization(self): "Ambiguous name 'test' matches 3 elements: ['test1', 'test2', 'test3']" ) + def test_object_matches_print_their_names(self): + "Regression: the message used to print raw objects instead of their names" + # WHEN matches are objects with unique names, as raised by EasyDynamicsList + matches = [ + Gaussian(name='SameName', unique_name='UniqueGaussian1'), + Gaussian(name='SameName', unique_name='UniqueGaussian2'), + ] + + # THEN + error = AmbiguousNameError('SameName', matches) + + # EXPECT the message names the matches instead of dumping object reprs + assert str(error) == ( + "Ambiguous name 'SameName' matches 2 elements: ['UniqueGaussian1', 'UniqueGaussian2']" + ) + def test_empty_matches(self): # WHEN name = 'unknown' diff --git a/tests/unit/easydynamics/test_import.py b/tests/unit/easydynamics/test_import.py deleted file mode 100644 index 11f87bdaf..000000000 --- a/tests/unit/easydynamics/test_import.py +++ /dev/null @@ -1,7 +0,0 @@ -# SPDX-FileCopyrightText: 2026 EasyScience contributors -# SPDX-License-Identifier: BSD-3-Clause - - -def test_import_easydynamics(): - # WHEN THEN EXPECT: importing raises no error - import easydynamics # ruff: ignore[unused-import] diff --git a/tests/unit/easydynamics/test_public_api.py b/tests/unit/easydynamics/test_public_api.py index fb81e81e4..e9dc71b4a 100644 --- a/tests/unit/easydynamics/test_public_api.py +++ b/tests/unit/easydynamics/test_public_api.py @@ -16,6 +16,7 @@ import pytest import easydynamics as edyn +from easydynamics.sample_model import Gaussian SUB_PACKAGES = [ 'easydynamics.analysis', @@ -23,6 +24,8 @@ 'easydynamics.convolution', 'easydynamics.experiment', 'easydynamics.sample_model', + 'easydynamics.sample_model.components', + 'easydynamics.sample_model.diffusion_model', 'easydynamics.settings', 'easydynamics.utils', ] @@ -31,6 +34,10 @@ class TestFrontDoor: + def test_import_easydynamics(self): + # WHEN THEN EXPECT: importing raises no error + import easydynamics # ruff: ignore[unused-import] + def test_everything_declared_is_importable(self): # THEN EXPECT no name in __all__ that cannot actually be reached missing = [name for name in edyn.__all__ if not hasattr(edyn, name)] @@ -47,9 +54,6 @@ def test_sub_package_exports_are_re_exported(self, module_name): assert missing == [], f'{module_name} exports not re-exported: {missing}' def test_re_exports_are_the_same_objects(self): - # WHEN - from easydynamics.sample_model import Gaussian - # THEN EXPECT the front door is an alias, not a copy assert edyn.Gaussian is Gaussian diff --git a/tests/unit/easydynamics/utils/test_detailed_balance.py b/tests/unit/easydynamics/utils/test_detailed_balance.py index 2d2284d36..f50d6ecad 100644 --- a/tests/unit/easydynamics/utils/test_detailed_balance.py +++ b/tests/unit/easydynamics/utils/test_detailed_balance.py @@ -17,29 +17,29 @@ class TestDetailedBalanceFactor: # Input validation tests def test_energy_unit_not_string_error(self): - # When + # WHEN energy = 2.0 T = 100 energy_unit = 5 - # Then Expect + # THEN EXPECT with pytest.raises(TypeError, match=r'energy_unit must be a string.'): detailed_balance_factor(energy, T, energy_unit=energy_unit) @pytest.mark.parametrize('temperature_unit', [5, 5.0, {}, []]) def test_temperature_unit_not_string_error(self, temperature_unit): - # When + # WHEN energy = 2.0 T = 100 - # Then Expect + # THEN EXPECT with pytest.raises(TypeError, match=r'temperature_unit must be a string.'): detailed_balance_factor(energy, T, temperature_unit=temperature_unit) def test_divide_by_temperature_not_bool_error(self): - # When + # WHEN energy = 2.0 T = 100 divide_by_temperature = 'yes' - # Then Expect + # THEN EXPECT with pytest.raises(TypeError, match=r'divide_by_temperature must be True or False.'): detailed_balance_factor(energy, T, divide_by_temperature=divide_by_temperature) @@ -61,11 +61,11 @@ def test_divide_by_temperature_not_bool_error(self): ], ) def test_energy_inputs(self, energy): - # When + # WHEN T = 100 - # Then + # THEN result = detailed_balance_factor(energy, T) - # Expect + # EXPECT if isinstance(energy, (np.ndarray)): energy_array = energy elif isinstance(energy, list): @@ -80,12 +80,12 @@ def test_energy_inputs(self, energy): np.testing.assert_allclose(result, expected, rtol=1e-5) def test_scipp_variable_input(self): - # When + # WHEN energy = sc.array(dims=['x'], values=[1.0, 2.0, 3.0], unit='meV') T = sc.scalar(value=100, unit='K') - # Then + # THEN result = detailed_balance_factor(energy, T) - # Expect + # EXPECT expected_values = ( np.array([1.0, 2.0, 3.0]) / (1 - np.exp(-np.array([1.0, 2.0, 3.0]) / (kB_meV_per_K * 100))) @@ -96,13 +96,58 @@ def test_scipp_variable_input(self): assert result.shape == (3,) np.testing.assert_allclose(result, expected_values, rtol=1e-5) + def test_dataarray_energy_input(self): + # WHEN + energy_values = np.array([1.0, 2.0, 3.0]) + data_array = sc.DataArray( + data=sc.array(dims=['energy'], values=np.zeros_like(energy_values)), + coords={ + 'energy': sc.array(dims=['energy'], values=energy_values, unit='meV'), + }, + ) + + # THEN + result = detailed_balance_factor(data_array, 100) + + # EXPECT the DataArray's single coordinate is used as the energy axis + expected = detailed_balance_factor(energy_values, 100) + np.testing.assert_allclose(result, expected) + + def test_dataarray_energy_with_multiple_coords_raises(self): + # WHEN + values = np.array([1.0, 2.0, 3.0]) + data_array = sc.DataArray( + data=sc.array(dims=['energy'], values=np.zeros_like(values)), + coords={ + 'energy': sc.array(dims=['energy'], values=values, unit='meV'), + 'other': sc.array(dims=['energy'], values=values, unit='meV'), + }, + ) + + # THEN EXPECT + with pytest.raises(ValueError, match='exactly one coordinate'): + detailed_balance_factor(data_array, 100) + + def test_two_dimensional_energy_raises(self): + # WHEN THEN EXPECT the documented ValueError, not a scipp DimensionError + with pytest.raises(ValueError, match='at most one-dimensional'): + detailed_balance_factor(np.ones((2, 2)), 100) + + def test_non_scalar_temperature_raises(self): + # WHEN + temperature = sc.array(dims=['temperature'], values=[100.0, 200.0], unit='K') + + # THEN EXPECT a clear error instead of a failure on `.value` + with pytest.raises(ValueError, match='temperature must be a single scalar value'): + detailed_balance_factor(np.array([1.0]), temperature) + def test_parameter_temperature(self): - # When + # WHEN energy = np.array([1.0, 2.0, 3.0]) T_param = Parameter(name='T', value=150, unit='K') - # Then + # THEN result = detailed_balance_factor(energy, T_param) - # Expect + # EXPECT expected = energy / (1 - np.exp(-energy / (kB_meV_per_K * 150))) / (kB_meV_per_K * 150) assert isinstance(result, np.ndarray) @@ -111,76 +156,76 @@ def test_parameter_temperature(self): # Physical edge cases def test_zero_temperature(self): - # When + # WHEN temperature = 0 energy = np.array([-1.0, 0.0, 1.0]) - # Then + # THEN result = detailed_balance_factor(energy, temperature, divide_by_temperature=False) - # Expect + # EXPECT expected = np.maximum(energy, 0.0) np.testing.assert_array_equal(result, expected) def test_zero_temperature_divide_by_T_error(self): - # When + # WHEN temperature = 0 energy = np.array([-1.0, 0.0, 1.0]) - # Then Expect + # THEN EXPECT with pytest.raises(ZeroDivisionError, match='Cannot divide by T when T = 0'): detailed_balance_factor(energy, temperature, divide_by_temperature=True) def test_zero_temperature_single_value(self): - # When + # WHEN temperature = 0 energy = 2.0 - # Then + # THEN result = detailed_balance_factor(energy, temperature, divide_by_temperature=False) - # Expect + # EXPECT expected = 2.0 assert result == expected def test_negative_temperature_raises(self): - # When Then Expect + # WHEN THEN EXPECT with pytest.raises(ValueError, match='Temperature must be non-negative'): detailed_balance_factor(1.0, -10) # Numerical tests def test_small_energy_limit(self): - # When + # WHEN T = 300 energy = np.array([1e-5, 1e-6, 1e-7, 1e-8, 1e-9]) - # Then + # THEN result = detailed_balance_factor(energy=energy, temperature=T, divide_by_temperature=False) - # Expect + # EXPECT x = energy / (kB_meV_per_K * T) expected = (1 + x / 2 + x**2 / 12) * (kB_meV_per_K * T) np.testing.assert_allclose(result, expected, rtol=1e-5) def test_large_energy_limit(self): - # When + # WHEN energy = np.linspace(1e2, 1e3, 5) T = 1 - # Then + # THEN result = detailed_balance_factor(energy=energy, temperature=T, divide_by_temperature=False) - # Expect + # EXPECT np.testing.assert_allclose(result, energy, atol=1e-10) def test_intermediate_energy(self): - # When + # WHEN energy = np.linspace(1, 10, 100) T = 100 - # Then + # THEN result = detailed_balance_factor(energy=energy, temperature=T, divide_by_temperature=False) - # Expect + # EXPECT expected = energy / (1 - np.exp(-energy / (kB_meV_per_K * T))) np.testing.assert_allclose(result, expected, rtol=1e-5) @pytest.mark.parametrize('divide_by_T', [True, False]) def test_detailed_balance_is_fulfilled(self, divide_by_T): # Detailed balance means DBF(E)/DBF(-E) = exp(E/(kB*T)) - # When + # WHEN T = 10 energy = np.linspace(0.01, 100, 101) - # Then + # THEN detailed_balance_positive = detailed_balance_factor( energy=energy, temperature=T, divide_by_temperature=divide_by_T ) @@ -189,7 +234,7 @@ def test_detailed_balance_is_fulfilled(self, divide_by_T): ) ratio = detailed_balance_positive / detailed_balance_negative - # Expect + # EXPECT expected_ratio = np.exp(energy / (kB_meV_per_K * T)) np.testing.assert_allclose(ratio, expected_ratio, rtol=1e-5) @@ -197,27 +242,27 @@ def test_detailed_balance_is_fulfilled(self, divide_by_T): 'energy_unit', ['microeV', sc.Unit('microeV')], ids=['str', 'scipp.Unit'] ) def test_energy_unit(self, energy_unit): - # When + # WHEN energy = np.linspace(1e3, 10 * 1e3, 100) T = 100 - # Then + # THEN result = detailed_balance_factor( energy=energy, temperature=T, divide_by_temperature=False, energy_unit=energy_unit, ) - # Expect + # EXPECT expected = energy / (1 - np.exp(-energy / 1000 / (kB_meV_per_K * T))) np.testing.assert_allclose(result, expected, rtol=1e-5) def test_energy_unit_warning(self): - # When + # WHEN energy = sc.linspace('energy', 1e3, 10 * 1e3, num=100, unit='microeV') energy_unit = 'meV' T = 100 - # Then + # THEN with pytest.warns( UserWarning, match='Input energy has unit [µμ]eV, but energy_unit was set to meV. Using [µμ]eV.', @@ -228,33 +273,33 @@ def test_energy_unit_warning(self): divide_by_temperature=False, energy_unit=energy_unit, ) - # Expect + # EXPECT expected = energy.values / (1 - np.exp(-energy.values / 1000 / (kB_meV_per_K * T))) np.testing.assert_allclose(result, expected, rtol=1e-5) @pytest.mark.parametrize('temperature_unit', ['mK', sc.Unit('mK')], ids=['str', 'scipp.Unit']) def test_temperature_unit(self, temperature_unit): - # When + # WHEN energy = np.linspace(1, 10, 100) temperature = 100 * 1000 temperature_unit = 'mK' - # Then + # THEN result = detailed_balance_factor( energy=energy, temperature=temperature, temperature_unit=temperature_unit, divide_by_temperature=False, ) - # Expect + # EXPECT expected = energy / (1 - np.exp(-energy / (kB_meV_per_K * temperature / 1000))) np.testing.assert_allclose(result, expected, rtol=1e-5) def test_temperature_unit_warning(self): - # When + # WHEN energy = np.linspace(1, 10, 100) temperature = sc.scalar(value=100, unit='mK') temperature_unit = 'K' - # Then + # THEN with pytest.warns( UserWarning, match='Input temperature has unit mK, but temperature_unit was set to K. Using mK.', @@ -265,18 +310,18 @@ def test_temperature_unit_warning(self): temperature_unit=temperature_unit, divide_by_temperature=False, ) - # Expect + # EXPECT expected = energy / (1 - np.exp(-energy / (kB_meV_per_K * 0.1))) np.testing.assert_allclose(result, expected, rtol=1e-5) def test_incompatible_energy_unit_raises(self): - # When + # WHEN energy = 2.0 T = 100 energy_unit = 'm' temperature_unit = 'K' - # Then Expect + # THEN EXPECT with pytest.raises( UnitError, match='The unit of energy is wrong', @@ -289,13 +334,13 @@ def test_incompatible_energy_unit_raises(self): ) def test_incompatible_temperature_unit_raises(self): - # When + # WHEN energy = 2.0 T = 100 energy_unit = 'meV' temperature_unit = 's' - # Then Expect + # THEN EXPECT with pytest.raises( UnitError, match='The unit of temperature is wrong', diff --git a/tests/unit/easydynamics/utils/test_fit_target.py b/tests/unit/easydynamics/utils/test_fit_target.py index bde22cb33..8c02658f7 100644 --- a/tests/unit/easydynamics/utils/test_fit_target.py +++ b/tests/unit/easydynamics/utils/test_fit_target.py @@ -8,51 +8,50 @@ from easydynamics.utils.fit_target import FitTarget -def test_fit_target_holds_prediction_metadata(): - # WHEN a FitTarget is created - target = FitTarget( - name='width', - dataset_key='Lorentzian width', - function=lambda x: x * 2, - label='DeltaLorentz width', - x_unit='1/angstrom', - y_unit='meV', - ) - # EXPECT its attributes to be preserved and the function callable - assert target.name == 'width' - assert target.dataset_key == 'Lorentzian width' - assert target.function(3) == 6 - assert target.label == 'DeltaLorentz width' - assert target.x_unit == '1/angstrom' - assert target.y_unit == 'meV' - - -def test_fit_target_allows_none_key_and_units(): - # WHEN a component-style FitTarget without a default key/units is created - target = FitTarget( - name='value', - dataset_key=None, - function=lambda x: x, - label='value', - x_unit=None, - y_unit=None, - ) - # EXPECT the optional fields to be None - assert target.dataset_key is None - assert target.x_unit is None - assert target.y_unit is None - - -def test_fit_target_is_frozen(): - # GIVEN a FitTarget - target = FitTarget( - name='value', - dataset_key=None, - function=lambda x: x, - label='value', - x_unit=None, - y_unit=None, - ) - # WHEN THEN EXPECT: it is immutable - with pytest.raises(FrozenInstanceError): - target.name = 'other' +class TestFitTarget: + def test_fit_target_holds_prediction_metadata(self): + # WHEN a FitTarget is created + target = FitTarget( + name='width', + dataset_key='Lorentzian width', + function=lambda x: x * 2, + label='DeltaLorentz width', + x_unit='1/angstrom', + y_unit='meV', + ) + # EXPECT its attributes to be preserved and the function callable + assert target.name == 'width' + assert target.dataset_key == 'Lorentzian width' + assert target.function(3) == 6 + assert target.label == 'DeltaLorentz width' + assert target.x_unit == '1/angstrom' + assert target.y_unit == 'meV' + + def test_fit_target_allows_none_key_and_units(self): + # WHEN a component-style FitTarget without a default key/units is created + target = FitTarget( + name='value', + dataset_key=None, + function=lambda x: x, + label='value', + x_unit=None, + y_unit=None, + ) + # EXPECT the optional fields to be None + assert target.dataset_key is None + assert target.x_unit is None + assert target.y_unit is None + + def test_fit_target_is_frozen(self): + # WHEN + target = FitTarget( + name='value', + dataset_key=None, + function=lambda x: x, + label='value', + x_unit=None, + y_unit=None, + ) + # THEN EXPECT: it is immutable + with pytest.raises(FrozenInstanceError): + target.name = 'other' diff --git a/tests/unit/easydynamics/utils/test_posterior_plotting.py b/tests/unit/easydynamics/utils/test_posterior_plotting.py index 1bb3bc4c3..892f1b6b2 100644 --- a/tests/unit/easydynamics/utils/test_posterior_plotting.py +++ b/tests/unit/easydynamics/utils/test_posterior_plotting.py @@ -162,6 +162,40 @@ def test_columns_share_limits_between_histogram_and_hexbin_panels(self, draws): column_limits = [grid[row, col].get_xlim() for row in range(col, 3)] assert all(limits == pytest.approx(column_limits[0]) for limits in column_limits) + ############# + # Scientific notation + ############# + + def test_shared_exponent_is_folded_into_the_label(self): + # WHEN the values are small enough that matplotlib factors out an exponent, which it parks + # on top of the axis label + draws = np.random.default_rng(0).normal(size=(200, 2)) * 1e-8 + 1.15e-8 + + # THEN + fig = plot_corner(draws=draws, names=['D', 'scale'], units=['m^2/s', '']) + + # EXPECT the exponent and the unit share one parenthetical, and the overlapping offset + # text is hidden + xlabel = fig.axes[-2].get_xlabel() + assert xlabel.startswith('D (1e') + assert 'm^2/s' in xlabel + assert not fig.axes[-2].xaxis.get_offset_text().get_visible() + + def test_shared_exponent_is_folded_into_the_y_label_too(self): + # WHEN the values are small enough that the left column's y axes also factor out an + # exponent + draws = np.random.default_rng(0).normal(size=(200, 2)) * 1e-8 + 1.15e-8 + + # THEN + fig = plot_corner(draws=draws, names=['D', 'scale'], units=['m^2/s', '']) + + # EXPECT the hexbin panel in the left column folds the exponent into its y label and + # hides the overlapping offset text + axis = fig.axes[2] + ylabel = axis.get_ylabel() + assert ylabel.startswith('scale (1e') + assert not axis.yaxis.get_offset_text().get_visible() + class TestPlotMarginal: @pytest.fixture @@ -578,20 +612,3 @@ def test_invalid_credible_interval_raises(self, arrays, interval): # THEN EXPECT with pytest.raises(ValueError, match='credible_interval'): predictive_with_slider(**arrays, credible_interval=interval) - - -class TestScientificNotation: - def test_shared_exponent_is_folded_into_the_label(self): - # WHEN the values are small enough that matplotlib factors out an exponent, which it parks - # on top of the axis label - draws = np.random.default_rng(0).normal(size=(200, 2)) * 1e-8 + 1.15e-8 - - # THEN - fig = plot_corner(draws=draws, names=['D', 'scale'], units=['m^2/s', '']) - - # EXPECT the exponent and the unit share one parenthetical, and the overlapping offset - # text is hidden - xlabel = fig.axes[-2].get_xlabel() - assert xlabel.startswith('D (1e') - assert 'm^2/s' in xlabel - assert not fig.axes[-2].xaxis.get_offset_text().get_visible() diff --git a/tests/unit/easydynamics/utils/test_utils.py b/tests/unit/easydynamics/utils/test_utils.py index a0bcaef97..023c4cfc5 100644 --- a/tests/unit/easydynamics/utils/test_utils.py +++ b/tests/unit/easydynamics/utils/test_utils.py @@ -54,6 +54,22 @@ def test_upper_bound_deferred_when_Q_is_none(self): # THEN EXPECT: a non-negative index is accepted; the bound check is deferred verify_Q_index(100, None) + @pytest.mark.parametrize('bool_index', [True, False], ids=['True', 'False']) + def test_bool_raises(self, bool_index): + # WHEN THEN EXPECT: bools are ints in Python, but Q_index=True must not mean index 1 + with pytest.raises(TypeError, match='Q_index must be an int'): + verify_Q_index(bool_index, None) + + def test_bool_raises_even_when_none_is_allowed(self): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match='Q_index must be an int or None'): + verify_Q_index(True, None, allow_none=True) + + def test_allow_none_rejects_non_int(self): + # WHEN THEN EXPECT: a non-int, non-None Q_index is rejected even when None is allowed + with pytest.raises(TypeError, match=r'Q_index must be an int or None'): + verify_Q_index('not an int', Q=None, allow_none=True) + class TestConvertValueUnit: def test_same_unit_returns_value_unchanged(self): @@ -302,40 +318,37 @@ def raise_import_error(*args, **kwargs): # ruff: ignore[unused-function-argumen assert _in_notebook() is False -def test_verify_Q_index_allow_none_rejects_non_int(): - # WHEN THEN EXPECT: a non-int, non-None Q_index is rejected even when None is allowed - with pytest.raises(TypeError, match=r'Q_index must be an int or None'): - verify_Q_index('not an int', Q=None, allow_none=True) - +class TestConvertParameterUnit: + def test_dependent_parameter_sets_desired_unit(self): + # WHEN converting the unit of a dependent parameter (cannot be converted directly) + param = Mock() + param.independent = False + convert_parameter_unit(param, 'meV') -def test_convert_parameter_unit_dependent_sets_desired_unit(): - # GIVEN a dependent parameter (cannot be converted directly) - param = Mock() - param.independent = False - # WHEN converting its unit - convert_parameter_unit(param, 'meV') - # EXPECT the desired unit is recorded instead of an in-place conversion - param.set_desired_unit.assert_called_once_with('meV') - param.convert_unit.assert_not_called() + # EXPECT the desired unit is recorded instead of an in-place conversion + param.set_desired_unit.assert_called_once_with('meV') + param.convert_unit.assert_not_called() -def test_energy_to_scipp_returns_variable_with_unit(): - # WHEN converting a numpy energy array - result = energy_to_scipp(np.array([1.0, 2.0, 3.0]), 'meV') - # EXPECT a scipp Variable on the 'energy' dimension with the given unit - assert isinstance(result, sc.Variable) - assert result.unit == sc.Unit('meV') - assert result.dims == ('energy',) - np.testing.assert_allclose(result.values, [1.0, 2.0, 3.0]) +class TestEnergyToScipp: + def test_returns_variable_with_unit(self): + # THEN + result = energy_to_scipp(np.array([1.0, 2.0, 3.0]), 'meV') + # EXPECT a scipp Variable on the 'energy' dimension with the given unit + assert isinstance(result, sc.Variable) + assert result.unit == sc.Unit('meV') + assert result.dims == ('energy',) + np.testing.assert_allclose(result.values, [1.0, 2.0, 3.0]) -def test_assert_valid_unit_rejects_non_unit_type(): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'unit must be a string or sc.Unit'): - _assert_valid_unit(123) +class TestAssertValidUnit: + def test_rejects_non_unit_type(self): + # THEN EXPECT + with pytest.raises(TypeError, match=r'unit must be a string or sc.Unit'): + _assert_valid_unit(123) -def test_assert_valid_unit_rejects_invalid_unit_string(): - # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'is not a valid scipp unit'): - _assert_valid_unit('not_a_real_unit') + def test_rejects_invalid_unit_string(self): + # THEN EXPECT + with pytest.raises(ValueError, match=r'is not a valid scipp unit'): + _assert_valid_unit('not_a_real_unit') diff --git a/tools/prefetch_tutorial_data.py b/tools/prefetch_tutorial_data.py index 839b1897a..12ebc8b84 100644 --- a/tools/prefetch_tutorial_data.py +++ b/tools/prefetch_tutorial_data.py @@ -23,8 +23,12 @@ TUTORIALS = Path(__file__).resolve().parent.parent / 'docs' / 'docs' / 'tutorials' -# Matches the pooch.retrieve(url=..., known_hash=...) calls the notebooks use, in either order. -URL_PATTERN = re.compile(r"url\s*=\s*f?['\"]([^'\"]+)['\"]") +# Matches one pooch.retrieve(...) call site; the notebooks keep these calls free of nested +# parentheses, so everything up to the first closing parenthesis is the argument list. +RETRIEVE_PATTERN = re.compile(r'pooch\.retrieve\s*\(([^)]*)\)') +# Matches the url=... and known_hash=... keyword arguments inside a single call. The optional +# ``f`` prefix on the URL is captured so templated URLs can be recognised and skipped. +URL_PATTERN = re.compile(r"url\s*=\s*(f?)['\"]([^'\"]+)['\"]") HASH_PATTERN = re.compile(r"known_hash\s*=\s*['\"]([^'\"]+)['\"]") @@ -32,26 +36,32 @@ def find_downloads() -> dict[str, str]: """ Collect the ``(url, known_hash)`` pairs the notebooks fetch. + Each ``pooch.retrieve(...)`` call site is parsed on its own, so a cell with several calls + cannot pair one call's URL with another call's hash. Calls whose URL is an f-string, or that + lack a literal ``known_hash``, are skipped rather than guessed at; the notebook will simply + fetch those itself. + Returns ------- dict[str, str] Mapping of URL to expected hash, deduplicated across notebooks. """ downloads: dict[str, str] = {} - for notebook in sorted(TUTORIALS.glob('*.ipynb')): + for notebook in sorted(TUTORIALS.rglob('*.ipynb')): + if '.ipynb_checkpoints' in notebook.parts: + continue cells = json.loads(notebook.read_text(encoding='utf-8'))['cells'] for cell in cells: if cell['cell_type'] != 'code': continue source = ''.join(cell['source']) - if 'pooch.retrieve' not in source: - continue - urls = URL_PATTERN.findall(source) - hashes = HASH_PATTERN.findall(source) - # Only pairs are usable; a templated URL without a literal hash is skipped rather than - # guessed at, and the notebook will simply fetch it itself. - for url, known_hash in zip(urls, hashes, strict=False): - downloads[url] = known_hash + for call in RETRIEVE_PATTERN.finditer(source): + arguments = call.group(1) + url_match = URL_PATTERN.search(arguments) + hash_match = HASH_PATTERN.search(arguments) + if url_match is None or hash_match is None or url_match.group(1) == 'f': + continue + downloads[url_match.group(2)] = hash_match.group(1) return downloads @@ -78,7 +88,7 @@ def main() -> int: name = url.rsplit('/', 1)[-1] try: pooch.retrieve(url=url, known_hash=known_hash) - except Exception as error: # noqa: BLE001 - report and continue, the notebook will retry + except Exception as error: # ruff: ignore[blind-except] - report and continue, the notebook will retry failures += 1 sys.stdout.write(f'could not prefetch {name}, leaving it to the notebook: {error}\n') else: