diff --git a/docs/docs/tutorials/components.ipynb b/docs/docs/tutorials/components.ipynb index c940c717..eafa8897 100644 --- a/docs/docs/tutorials/components.ipynb +++ b/docs/docs/tutorials/components.ipynb @@ -151,6 +151,15 @@ "expr = sm.ExpressionComponent(\n", " 'A * exp(-(x - x0)**2 / (2*sigma**2)) +B*sin(2*pi*x/period)',\n", " parameters={'A': 10, 'x0': 0, 'sigma': 1},\n", + " parameter_units={\n", + " 'A': 'dimensionless',\n", + " 'x0': 'meV',\n", + " 'sigma': 'meV',\n", + " 'B': 'dimensionless',\n", + " 'period': 'meV',\n", + " },\n", + " x_unit='meV',\n", + " y_unit='dimensionless',\n", ")\n", "\n", "expr.A = 5\n", diff --git a/src/easydynamics/sample_model/components/expression_component.py b/src/easydynamics/sample_model/components/expression_component.py index eef67ff6..8db32c4d 100644 --- a/src/easydynamics/sample_model/components/expression_component.py +++ b/src/easydynamics/sample_model/components/expression_component.py @@ -4,19 +4,24 @@ from __future__ import annotations import warnings +from copy import copy from typing import TYPE_CHECKING from typing import ClassVar +import scipp as sc import sympy as sp +from easyscience.variable import DescriptorNumber from easyscience.variable import Parameter from scipy.special import erf if TYPE_CHECKING: import numpy as np - import scipp as sc from easydynamics.sample_model.components.model_component import ModelComponent from easydynamics.utils.utils import Numeric +from easydynamics.utils.utils import _assert_valid_unit +from easydynamics.utils.utils import hbar +from easydynamics.utils.utils import kb class ExpressionComponent(ModelComponent): @@ -55,6 +60,36 @@ class ExpressionComponent(ModelComponent): expr.A = 5 expr.sigma = 0.5 ``` + + **Giving parameters units** + + 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: + ```python + expr = sm.ExpressionComponent( + 'A * exp(-(x - x0)**2 / (2*sigma**2))', + parameters={'A': 10, 'x0': 0, 'sigma': 1}, + parameter_units={'A': '1/meV', 'x0': 'meV', 'sigma': 'meV'}, + y_unit='1/meV', + ) + expr.set_unit('A', '1/meV') + ``` + + **Physical constants** + + The symbols ``hbar`` (in meV*s) and ``kb`` (in meV/K) are provided automatically as read-only + constants (DescriptorNumbers) when they appear in the expression: + ```python + boltzmann = sm.ExpressionComponent( + 'exp(-x / (kb * T))', + parameters={'T': 300.0}, + parameter_units={'T': 'K'}, + ) + ``` + Use e.g. ``boltzmann.kb.convert_unit('eV/K')`` to work in another unit (this rescales the + value, unlike ``set_unit``). """ _ALLOWED_FUNCS: ClassVar[dict[str, object]] = { @@ -95,10 +130,28 @@ class ExpressionComponent(ModelComponent): _RESERVED_NAMES: ClassVar[dict[str, object]] = {'x'} + # Physical constants provided automatically as read-only DescriptorNumbers when their + # symbol appears in the expression: name -> (source constant from utils, default unit). + _PHYSICAL_CONSTANTS: ClassVar[dict[str, tuple[DescriptorNumber, str]]] = { + 'hbar': (hbar, 'meV*s'), + 'kb': (kb, 'meV/K'), + } + + # Sympy functions that preserve the unit of their argument; all other allowed functions + # require a dimensionless argument and return a dimensionless result. + _UNIT_PRESERVING_FUNCS: ClassVar[tuple] = (sp.Abs, sp.floor, sp.ceiling) + # Functions that accept any unit but return a dimensionless result. + _DIMENSIONLESS_RESULT_FUNCS: ClassVar[tuple] = (sp.sign,) + + # parameter_units is applied at construction only; the resulting units are captured in + # the serialized Parameter dicts, so the argument is skipped during serialization. + _REDIRECT: ClassVar[dict[str, object]] = {'parameter_units': None} + def __init__( self, expression: str, parameters: dict[str, Numeric] | None = None, + parameter_units: dict[str, str | sc.Unit] | None = None, x_unit: str | sc.Unit = 'meV', y_unit: str | sc.Unit = 'dimensionless', name: str = 'Expression', @@ -111,9 +164,17 @@ def __init__( Parameters ---------- expression : str - The symbolic expression as a string. Must contain 'x' as the independent variable. + 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*. parameters : dict[str, Numeric] | None, default=None - Dictionary of parameter names and their initial values. + 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. x_unit : str | sc.Unit, default='meV' Unit of the x-axis. y_unit : str | sc.Unit, default='dimensionless' @@ -128,9 +189,10 @@ def __init__( Raises ------ ValueError - If the expression is invalid or does not contain 'x'. + If the expression is invalid or does not contain 'x', or if parameter_units names a + parameter that is not in the expression. TypeError - If any parameter value is not numeric. + If any parameter value is not numeric, or if parameter_units is not a dictionary. """ super().__init__( x_unit=x_unit, @@ -198,12 +260,19 @@ def __init__( ) parameters = parameters or {} self._parameters: dict[str, Parameter] = {} + self._constants: dict[str, DescriptorNumber] = {} self._symbol_names = symbol_names for name in self._symbol_names: if name in self._RESERVED_NAMES: continue + # 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: + self._constants[name] = self._create_physical_constant(name) + continue + value = parameters.get(name, 1.0) if isinstance(value, Parameter): self._parameters[name] = value @@ -213,8 +282,32 @@ def __init__( self._parameters[name] = Parameter( name=name, value=value, - unit=self.x_unit, + unit='dimensionless', + ) + + if parameter_units is not None: + if not isinstance(parameter_units, dict): + raise TypeError( + f'parameter_units must be None or a dictionary, ' + f'got {type(parameter_units).__name__}' ) + constant_names = sorted(set(parameter_units) & set(self._constants)) + if constant_names: + raise ValueError( + f'parameter_units given for physical constant(s): ' + f'{", ".join(constant_names)}. Use e.g. ' + f'expr.{constant_names[0]}.convert_unit(...) to change the unit of a ' + f'constant (this rescales its value).' + ) + unknown_names = sorted(set(parameter_units) - set(self._parameters)) + if unknown_names: + raise ValueError( + f'parameter_units given for unknown parameter(s): {", ".join(unknown_names)}' + ) + for parameter_name, parameter_unit in parameter_units.items(): + # Relabel without the per-call output-unit warning; the consistency check + # runs once below, after all units are applied. + self._relabel_parameter_unit(parameter_name, parameter_unit) # Create numerical function ordered_symbols = [sp.Symbol(name) for name in self._symbol_names] @@ -224,6 +317,8 @@ def __init__( modules=[{'erf': erf}, 'numpy'], ) + self._warn_if_output_unit_mismatch() + # ------------------------- # Properties # ------------------------- @@ -263,8 +358,8 @@ def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndar Evaluate the expression for the given x values. Unit conversion of parameters is not supported for ExpressionComponent. If x_vals is - expressed in a different unit than x_unit, a warning is issued and the values are used - as-is. + expressed in a different unit than x_unit, a UnitError is raised — silently evaluating the + expression with wrongly-scaled x would corrupt the result. Parameters ---------- @@ -273,24 +368,33 @@ def _evaluate_values(self, x_vals: np.ndarray, eval_unit: str | None) -> np.ndar eval_unit : str | None The unit of x_vals. + Raises + ------ + sc.UnitError + If eval_unit differs from the component's x_unit. + Returns ------- np.ndarray Evaluated results. """ - if eval_unit is not None and self.x_unit is not None and eval_unit != self.x_unit: - warnings.warn( + if ( + eval_unit is not None + and self.x_unit is not None + and sc.Unit(eval_unit) != sc.Unit(self.x_unit) + ): + raise sc.UnitError( f'Input x has unit {eval_unit} but {self.__class__.__name__} has ' - f'x_unit {self.x_unit}. ExpressionComponent cannot auto-convert parameters. ' - 'x values are used as-is.', - UserWarning, - stacklevel=3, + f'x_unit {self.x_unit}. ExpressionComponent cannot auto-convert its parameters; ' + f'convert x to {self.x_unit} before evaluating.' ) args = [] for name in self._symbol_names: if name == 'x': args.append(x_vals) + elif name in self._constants: + args.append(self._constants[name].value) else: args.append(self._parameters[name].value) @@ -307,11 +411,299 @@ def get_all_variables(self) -> list[Parameter]: """ return list(self._parameters.values()) + def set_unit(self, name: str, unit: str | sc.Unit) -> None: + """ + Set the unit of a parameter without rescaling its value. + + 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. + + Parameters + ---------- + name : str + Name of the parameter whose unit to set. + unit : str | sc.Unit + The new unit. + """ + self._relabel_parameter_unit(name, unit) + self._warn_if_output_unit_mismatch() + + def _relabel_parameter_unit(self, name: str, unit: str | sc.Unit) -> None: + """ + Validate and apply a unit relabel to a parameter (see set_unit). + + Unit validation is delegated to ``_assert_valid_unit``, which raises ``ValueError`` for a + string that is not a valid scipp unit. + + Parameters + ---------- + name : str + Name of the parameter whose unit to set. + unit : str | sc.Unit + The new unit. + + Raises + ------ + TypeError + If name is not a string or unit is not a string or sc.Unit. + KeyError + If no parameter with the given name exists. + AttributeError + If the parameter is a physical constant or a dependent parameter. + """ + if not isinstance(name, str): + raise TypeError(f'name must be a string, got {type(name).__name__}') + if '_constants' in self.__dict__ and name in self._constants: + raise AttributeError( + f"'{name}' is a physical constant. Use expr.{name}.convert_unit(...) to change " + f'its unit (this rescales its value).' + ) + if '_parameters' not in self.__dict__ or name not in self._parameters: + raise KeyError(f"No parameter named '{name}' in this {self.__class__.__name__}.") + _assert_valid_unit(unit) + new_unit = sc.Unit(str(unit)) + + param = self._parameters[name] + if not param.independent: + raise AttributeError( + f"Cannot set the unit of dependent parameter '{name}'; its unit is controlled " + f'by the dependency expression.' + ) + + # Parameter.unit is read-only and convert_unit rescales the value, so a pure relabel + # has to swap the underlying scipp scalars (value and bounds) directly. + param._scalar = sc.scalar( # noqa: SLF001 + param.value, + unit=new_unit, + variance=param._scalar.variance, # noqa: SLF001 + ) + param._min = sc.scalar(param._min.value, unit=new_unit) # noqa: SLF001 + param._max = sc.scalar(param._max.value, unit=new_unit) # noqa: SLF001 + + @classmethod + def _create_physical_constant(cls, name: str) -> DescriptorNumber: + """ + Create a DescriptorNumber for a supported physical constant. + + The constant is a per-instance copy of the shared constant in ``easydynamics.utils``, + converted to its default unit, so converting its unit on one component does not affect + others. + + Parameters + ---------- + name : str + The constant's symbol name (a key of ``_PHYSICAL_CONSTANTS``). + + Returns + ------- + DescriptorNumber + The constant's value in its default unit. + """ + source, default_unit = cls._PHYSICAL_CONSTANTS[name] + constant = copy(source) + constant.convert_unit(default_unit) + return constant + + @property + def constants(self) -> dict[str, DescriptorNumber]: + """ + Get the physical constants used by the expression. + + Returns + ------- + dict[str, DescriptorNumber] + The automatically provided constants (e.g. ``hbar``, ``kb``) keyed by symbol name. + """ + return dict(self._constants) + + @property + def output_unit(self) -> str: + """ + Get the unit of the evaluated expression, derived from x_unit and the parameter units. + + The unit is propagated through the expression tree: addition requires compatible units, + multiplication and powers combine units, and functions like ``exp`` or ``sin`` require a + dimensionless argument. Propagation raises ``sc.UnitError`` if the expression is not + unit-consistent (e.g. adding meV to a dimensionless quantity, or taking ``exp`` of a + quantity with a unit). + + Returns + ------- + str + The unit of the evaluated expression. + """ + return str(self._canonical_unit(self._propagate_unit(self._expr))) + + def _canonical_unit(self, unit: sc.Unit) -> sc.Unit: + """ + Replace a unit by an equal, cleaner-printing candidate where possible. + + Unit algebra accumulates tiny floating-point scale factors (e.g. ``meV * 1/meV`` prints as + ``0.999999999999999889`` rather than ``dimensionless``); scipp's unit equality still + matches, so map such units onto the common candidates for readable output. + + Parameters + ---------- + unit : sc.Unit + The unit to canonicalize. + + Returns + ------- + sc.Unit + An equal unit with a cleaner string representation, or the input unit unchanged. + """ + candidate_strings = ['dimensionless', self.y_unit, self.x_unit] + for candidate_string in candidate_strings: + if candidate_string is None: + continue + candidate = sc.Unit(candidate_string) + if unit == candidate: + return candidate + return unit + + def _propagate_unit(self, node: sp.Basic) -> sc.Unit: + """ + Recursively derive the unit of a sympy expression node. + + Parameters + ---------- + node : sp.Basic + The sympy expression node to derive the unit of. + + Raises + ------ + sc.UnitError + If the node mixes incompatible units, applies a function that requires a dimensionless + argument to a quantity with a unit, or raises a quantity with a unit to a symbolic or + unsupported power. + + Returns + ------- + sc.Unit + The unit of the node. + """ + dimensionless = sc.Unit('dimensionless') + + if isinstance(node, sp.Symbol): + name = str(node) + if name == 'x': + return sc.Unit(self.x_unit) if self.x_unit is not None else dimensionless + if name in self._constants: + return sc.Unit(str(self._constants[name].unit)) + return sc.Unit(str(self._parameters[name].unit)) + + if node.is_number: + return dimensionless + + if isinstance(node, sp.Add): + units = [self._propagate_unit(argument) for argument in node.args] + for unit in units[1:]: + try: + sc.to_unit(sc.scalar(1.0, unit=unit), units[0]) + except sc.UnitError as e: + raise sc.UnitError( + f'The expression adds quantities with incompatible units ' + f'{units[0]} and {unit}.' + ) from e + return units[0] + + if isinstance(node, sp.Mul): + result = dimensionless + for argument in node.args: + result = result * self._propagate_unit(argument) + return result + + if isinstance(node, sp.Pow): + base_unit = self._propagate_unit(node.base) + exponent_unit = self._propagate_unit(node.exp) + if exponent_unit != dimensionless: + raise sc.UnitError(f'An exponent must be dimensionless, got {exponent_unit}.') + if base_unit == dimensionless: + return dimensionless + if not node.exp.is_number: + raise sc.UnitError( + f'Cannot determine units: quantity with unit {base_unit} is raised to a ' + f'symbolic power.' + ) + exponent = float(node.exp) + if exponent == int(exponent): + return base_unit ** int(exponent) + if 2 * exponent == int(2 * exponent): + # Half-integer power: the square root of an integer power. + return sc.sqrt(sc.scalar(1.0, unit=base_unit ** int(2 * exponent))).unit + raise sc.UnitError( + f'Cannot determine units: quantity with unit {base_unit} is raised to the ' + f'power {exponent}.' + ) + + if isinstance(node, self._UNIT_PRESERVING_FUNCS): + return self._propagate_unit(node.args[0]) + + if isinstance(node, self._DIMENSIONLESS_RESULT_FUNCS): + self._propagate_unit(node.args[0]) + return dimensionless + + if isinstance(node, sp.Function): + for argument in node.args: + argument_unit = self._propagate_unit(argument) + try: + sc.to_unit(sc.scalar(1.0, unit=argument_unit), dimensionless) + except sc.UnitError as e: + raise sc.UnitError( + f'{node.func.__name__} requires a dimensionless argument, ' + f'got {argument_unit}.' + ) from e + return dimensionless + + raise sc.UnitError( + f'Cannot determine units for expression node {node} of type {type(node).__name__}.' + ) + + def _warn_if_output_unit_mismatch(self) -> None: + """ + Warn if the expression's output unit does not match 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. + """ + units_in_use = bool(self._constants) or any( + str(parameter.unit) != 'dimensionless' for parameter in self._parameters.values() + ) + if not units_in_use: + return + + try: + output_unit = sc.Unit(self.output_unit) + except sc.UnitError as e: + warnings.warn( + f'The expression is not unit-consistent: {e}', + UserWarning, + stacklevel=3, + ) + return + + y_unit = sc.Unit(self.y_unit) if self.y_unit is not None else sc.Unit('dimensionless') + if output_unit != y_unit: + 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.', + UserWarning, + stacklevel=3, + ) + def convert_x_unit(self, _new_unit: str | sc.Unit) -> None: """ Convert the x-axis unit of the expression. - Unit conversion is not implemented for ExpressionComponent. + Unit conversion is not implemented for ExpressionComponent. Should it ever be needed, the + viable path is dimensional analysis on the parameter units: for each parameter, determine + the power n of the x-dimension in its unit and rescale its value by the x-unit conversion + factor to the power n (the generalization of Polynomial's power-law rescaling). This only + works when x_unit has a single unambiguous dimension. Parameters ---------- @@ -323,14 +715,14 @@ def convert_x_unit(self, _new_unit: str | sc.Unit) -> None: NotImplementedError Always raised to indicate unit conversion is not supported. """ - # TODO: Generalize x_unit conversion for ExpressionComponent. Not tackled in this PR. # noqa: FIX002 TD002 TD003 raise NotImplementedError('Unit conversion is not implemented for ExpressionComponent') def convert_y_unit(self, _new_unit: str | sc.Unit) -> None: """ Convert the y-axis unit of the expression. - Unit conversion is not implemented for ExpressionComponent. + Unit conversion is not implemented for ExpressionComponent. See convert_x_unit for the + approach that would make it possible. Parameters ---------- @@ -342,21 +734,20 @@ def convert_y_unit(self, _new_unit: str | sc.Unit) -> None: NotImplementedError Always raised to indicate unit conversion is not supported. """ - # TODO: Generalize y_unit conversion for ExpressionComponent. Not tackled in this PR. # noqa: FIX002 TD002 TD003 raise NotImplementedError('Unit conversion is not implemented for ExpressionComponent') # ------------------------- # dunder methods # ------------------------- - def __getattr__(self, name: str) -> Parameter: + def __getattr__(self, name: str) -> Parameter | DescriptorNumber: """ - Allow access to parameters as attributes. + Allow access to parameters and physical constants as attributes. Parameters ---------- name : str - Name of the parameter to access. + Name of the parameter or constant to access. Raises ------ @@ -365,11 +756,13 @@ def __getattr__(self, name: str) -> Parameter: Returns ------- - Parameter - The parameter with the given name. + Parameter | DescriptorNumber + The parameter or constant with the given name. """ if '_parameters' in self.__dict__ and name in self._parameters: return self._parameters[name] + if '_constants' in self.__dict__ and name in self._constants: + return self._constants[name] raise AttributeError(f"{self.__class__.__name__} has no attribute '{name}'") def __setattr__(self, name: str, value: Numeric) -> None: @@ -385,9 +778,13 @@ def __setattr__(self, name: str, value: Numeric) -> None: Raises ------ + AttributeError + If the name refers to a physical constant. TypeError If the value is not numeric. """ + if '_constants' in self.__dict__ and name in self._constants: + raise AttributeError(f"'{name}' is a physical constant and cannot be set.") if '_parameters' in self.__dict__ and name in self._parameters: param = self._parameters[name] if not isinstance(value, Numeric): @@ -398,14 +795,14 @@ def __setattr__(self, name: str, value: Numeric) -> None: def __dir__(self) -> list[str]: """ - Include parameter names in dir() output for better IDE support. + Include parameter and constant names in dir() output for better IDE support. Returns ------- list[str] - List of attribute names, including parameters. + List of attribute names, including parameters and constants. """ - return super().__dir__() + list(self._parameters.keys()) + return super().__dir__() + list(self._parameters.keys()) + list(self._constants.keys()) def __repr__(self) -> str: """ @@ -417,9 +814,15 @@ def __repr__(self) -> str: String representation of the ExpressionComponent. """ param_str = ', '.join(f'{k}={v.value}' for k, v in self._parameters.items()) + constants_str = '' + if self._constants: + constants_str = ',\n constants={ ' + ', '.join( + f'{k}={v.value} {v.unit}' for k, v in self._constants.items() + ) + constants_str += ' }' return ( f'{self.__class__.__name__}(name={self.name}, display_name={self.display_name}, ' f'x_unit={self.x_unit}, y_unit={self.y_unit},\n' f" expr='{self._expression_str}',\n" - f' parameters={{ {param_str} }} )' + f' parameters={{ {param_str} }}{constants_str} )' ) 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 6cf81147..ba34ff5f 100644 --- a/tests/unit/easydynamics/sample_model/components/test_expression_component.py +++ b/tests/unit/easydynamics/sample_model/components/test_expression_component.py @@ -1,14 +1,18 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import warnings from copy import copy import numpy as np import pytest import scipp as sc +from easyscience.variable import DescriptorNumber from easyscience.variable import Parameter from easydynamics.sample_model import ExpressionComponent +from easydynamics.sample_model import Gaussian +from easydynamics.sample_model import Lorentzian class TestExpressionComponent: @@ -45,6 +49,50 @@ def test_init_with_parameter(self): # EXPECT assert expr.A.value == pytest.approx(3.0) + def test_parameters_are_dimensionless_by_default(self, expr: ExpressionComponent): + # WHEN THEN EXPECT + assert str(expr.A.unit) == 'dimensionless' + assert str(expr.x0.unit) == 'dimensionless' + assert str(expr.sigma.unit) == 'dimensionless' + + def test_init_with_parameter_units(self): + # WHEN THEN + expr = ExpressionComponent( + 'A * exp(-(x - x0)**2 / (2*sigma**2))', + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units={'x0': 'meV', 'sigma': 'meV'}, + ) + + # EXPECT: the named parameters get their units, values unchanged; others stay + # dimensionless + assert str(expr.x0.unit) == 'meV' + assert str(expr.sigma.unit) == 'meV' + assert str(expr.A.unit) == 'dimensionless' + assert expr.x0.value == pytest.approx(0.5) + assert expr.sigma.value == pytest.approx(0.6) + + def test_init_parameter_units_overrides_parameter_instance_unit(self): + # WHEN: a Parameter instance in eV, but parameter_units requests meV + A = Parameter('A', 3.0, unit='eV') + + # THEN: A*x evaluates to meV**2, which mismatches the default dimensionless y_unit + with pytest.warns(UserWarning, match='does not match'): + expr = ExpressionComponent('A * x', parameters={'A': A}, parameter_units={'A': 'meV'}) + + # EXPECT: the unit is relabelled without rescaling the value + assert str(expr.A.unit) == 'meV' + assert expr.A.value == pytest.approx(3.0) + + def test_init_parameter_units_unknown_name_raises(self): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='unknown parameter'): + ExpressionComponent('A * x', parameter_units={'B': 'meV'}) + + def test_init_parameter_units_not_a_dict_raises(self): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match='parameter_units must be None or a dictionary'): + ExpressionComponent('A * x', parameter_units='meV') + def test_invalid_expression_raises(self): # WHEN THEN EXPECT with pytest.raises(ValueError, match='Invalid expression'): @@ -132,6 +180,90 @@ def test_expression_is_read_only(self, expr: ExpressionComponent): with pytest.raises(AttributeError, match='cannot be changed'): expr.expression = 'x' + def test_set_unit_relabels_without_rescaling(self, expr: ExpressionComponent): + # WHEN + expr.sigma.min = 0.1 + expr.sigma.max = 10.0 + + # THEN: relabelling only sigma leaves the expression unit-inconsistent, which warns + with pytest.warns(UserWarning, match='not unit-consistent'): + expr.set_unit('sigma', 'meV') + + # EXPECT: unit relabelled; value and bounds keep their numeric values + assert str(expr.sigma.unit) == 'meV' + assert expr.sigma.value == pytest.approx(0.6) + assert expr.sigma.min == pytest.approx(0.1) + assert expr.sigma.max == pytest.approx(10.0) + + def test_set_unit_accepts_scipp_unit(self, expr: ExpressionComponent): + # WHEN THEN: relabelling only x0 leaves the expression unit-inconsistent, which warns + with pytest.warns(UserWarning, match='not unit-consistent'): + expr.set_unit('x0', sc.Unit('meV')) + + # EXPECT + assert str(expr.x0.unit) == 'meV' + + def test_set_unit_leaves_other_parameters_untouched(self, expr: ExpressionComponent): + # WHEN THEN + with pytest.warns(UserWarning, match='not unit-consistent'): + expr.set_unit('x0', 'meV') + + # EXPECT + assert str(expr.A.unit) == 'dimensionless' + assert str(expr.sigma.unit) == 'dimensionless' + + def test_set_unit_unknown_parameter_raises(self, expr: ExpressionComponent): + # WHEN THEN EXPECT + with pytest.raises(KeyError, match="No parameter named 'B'"): + expr.set_unit('B', 'meV') + + def test_set_unit_invalid_unit_raises(self, expr: ExpressionComponent): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='not a valid scipp unit'): + expr.set_unit('A', 'not_a_unit') + + @pytest.mark.parametrize( + 'name, unit, expected_message', + [ + (123, 'meV', 'name must be a string'), + ('A', 123, 'unit must be a string or sc.Unit'), + ], + ids=['non_string_name', 'non_unit_unit'], + ) + def test_set_unit_type_validation( + self, expr: ExpressionComponent, name, unit, expected_message + ): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match=expected_message): + expr.set_unit(name, unit) + + def test_set_unit_dependent_parameter_raises(self, expr: ExpressionComponent): + # WHEN: make A dependent on sigma + expr.A.make_dependent_on( + dependency_expression='2 * sigma', + dependency_map={'sigma': expr.sigma}, + ) + + # THEN EXPECT + with pytest.raises(AttributeError, match='dependent parameter'): + expr.set_unit('A', 'meV') + + def test_parameter_units_survive_copy(self): + # WHEN + expr = ExpressionComponent( + 'A * x + B', + parameters={'A': 2.0, 'B': 3.0}, + parameter_units={'A': '1/meV'}, + ) + + # THEN + expr_copy = copy(expr) + + # EXPECT: per-parameter units round-trip through to_dict/from_dict + assert str(expr_copy.A.unit) == '1/meV' + assert str(expr_copy.B.unit) == 'dimensionless' + assert expr_copy.A.value == pytest.approx(2.0) + def test_convert_x_unit_not_implemented(self, expr: ExpressionComponent): # WHEN THEN EXPECT with pytest.raises(NotImplementedError, match='not implemented'): @@ -228,10 +360,423 @@ def test_erf(self): np.testing.assert_allclose(result, expected, rtol=1e-5) -def test_evaluate_warns_when_input_unit_differs_from_x_unit(): +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 warning - with pytest.warns(UserWarning, match=r'cannot auto-convert parameters'): + # 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.""" + + @pytest.fixture + def gaussian_expr(self): + return ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units=GAUSSIAN_UNITS, + x_unit='meV', + ) + + def test_matches_gaussian_component(self, gaussian_expr): + # WHEN + gaussian = Gaussian(area=2.0, center=0.5, width=0.6, x_unit='meV') + x = np.linspace(-3, 3, 101) + + # THEN EXPECT: same values on a plain numpy grid + np.testing.assert_allclose(gaussian_expr.evaluate(x), gaussian.evaluate(x)) + + def test_matches_gaussian_component_scipp_input_and_output(self, gaussian_expr): + # WHEN + gaussian = Gaussian(area=2.0, center=0.5, width=0.6, x_unit='meV') + x = sc.linspace('energy', -3.0, 3.0, 101, unit='meV') + + # THEN + expr_result = gaussian_expr.evaluate(x, output='scipp') + gaussian_result = gaussian.evaluate(x, output='scipp') + + # EXPECT: same values and the same output unit + np.testing.assert_allclose(expr_result.values, gaussian_result.values) + assert expr_result.unit == gaussian_result.unit + + def test_matches_gaussian_component_in_other_unit(self): + # WHEN: the built-in Gaussian auto-converts its meV parameters to a ueV grid; the + # expression is given the physically identical parameters expressed in ueV directly. + gaussian = Gaussian(area=2.0, center=0.5, width=0.6, x_unit='meV') + expr = ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2000.0, 'x0': 500.0, 'sigma': 600.0}, + parameter_units={'A': 'ueV', 'x0': 'ueV', 'sigma': 'ueV'}, + x_unit='ueV', + ) + x = sc.linspace('energy', -3000.0, 3000.0, 101, unit='ueV') + + # THEN EXPECT: identical physical results (the Gaussian output is dimensionless, so no + # y-scale factor is involved) + np.testing.assert_allclose(expr.evaluate(x), gaussian.evaluate(x)) + + def test_matches_lorentzian_component(self): + # WHEN + lorentzian = Lorentzian(area=2.0, center=0.5, width=0.3, x_unit='meV') + # Note: 'gamma' would collide with sympy's gamma function, so use 'hwhm' + expr = ExpressionComponent( + 'A / pi * hwhm / ((x - x0)**2 + hwhm**2)', + parameters={'A': 2.0, 'x0': 0.5, 'hwhm': 0.3}, + parameter_units={'A': 'meV', 'x0': 'meV', 'hwhm': 'meV'}, + x_unit='meV', + ) + x = np.linspace(-3, 3, 101) + + # THEN EXPECT + np.testing.assert_allclose(expr.evaluate(x), lorentzian.evaluate(x)) + + def test_consistent_units_do_not_warn(self): + # WHEN THEN EXPECT: a fully unit-consistent expression constructs without warnings + with warnings.catch_warnings(): + warnings.simplefilter('error') + ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units=GAUSSIAN_UNITS, + x_unit='meV', + ) + + def test_unit_agnostic_expression_does_not_warn(self): + # WHEN THEN EXPECT: without any parameter units the consistency check stays silent, + # even though x_unit is meV and the parameters are dimensionless + with warnings.catch_warnings(): + warnings.simplefilter('error') + ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + x_unit='meV', + ) + + +class TestExpressionComponentOutputUnit: + def test_output_unit_gaussian_is_dimensionless(self): + # WHEN: area in meV divided by sigma in meV + expr = ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units=GAUSSIAN_UNITS, + x_unit='meV', + ) + + # THEN EXPECT + assert expr.output_unit == 'dimensionless' + + def test_output_unit_with_dimensionless_area_is_one_over_x(self): + # WHEN: dimensionless amplitude over a meV width gives 1/meV + with pytest.warns(UserWarning, match='does not match'): + expr = ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units={'x0': 'meV', 'sigma': 'meV'}, + x_unit='meV', + ) + + # THEN EXPECT + assert sc.Unit(expr.output_unit) == sc.Unit('1/meV') + + def test_output_unit_matching_y_unit_does_not_warn(self): + # WHEN THEN EXPECT: declaring the matching y_unit silences the mismatch warning + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent( + GAUSSIAN_EXPRESSION, + parameters={'A': 2.0, 'x0': 0.5, 'sigma': 0.6}, + parameter_units={'x0': 'meV', 'sigma': 'meV'}, + x_unit='meV', + y_unit='1/meV', + ) + assert sc.Unit(expr.output_unit) == sc.Unit(expr.y_unit) + + def test_output_unit_preserved_by_abs(self): + # WHEN + with pytest.warns(UserWarning, match='does not match'): + expr = ExpressionComponent( + 'abs(x - x0)', parameters={'x0': 0.5}, parameter_units={'x0': 'meV'} + ) + + # THEN EXPECT + assert sc.Unit(expr.output_unit) == sc.Unit('meV') + + def test_output_unit_sign_is_dimensionless(self): + # WHEN + expr = ExpressionComponent( + 'A * sign(x - x0)', parameters={'A': 1.0, 'x0': 0.5}, parameter_units={'x0': 'meV'} + ) + + # THEN EXPECT + assert expr.output_unit == 'dimensionless' + + def test_output_unit_sqrt_of_squared_quantity(self): + # WHEN: sqrt(sigma**2) is a half-integer power of a quantity with units + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent( + 'x / sqrt(2*pi*sigma**2) + offset', + parameters={'sigma': 0.5, 'offset': 1.0}, + parameter_units={'sigma': 'meV'}, + x_unit='meV', + ) + + # THEN EXPECT: meV / sqrt(meV**2) + dimensionless = dimensionless + assert expr.output_unit == 'dimensionless' + + def test_output_unit_incompatible_addition_raises(self): + # WHEN + with pytest.warns(UserWarning, match='not unit-consistent'): + expr = ExpressionComponent( + 'x + A', parameters={'A': 1.0}, parameter_units={'A': 'K'}, x_unit='meV' + ) + + # THEN EXPECT + with pytest.raises(sc.UnitError, match='incompatible units'): + _ = expr.output_unit + + def test_output_unit_exp_of_dimensional_quantity_raises(self): + # WHEN + with pytest.warns(UserWarning, match='not unit-consistent'): + expr = ExpressionComponent( + 'exp(x * A)', parameters={'A': 1.0}, parameter_units={'A': 'meV'}, x_unit='meV' + ) + + # THEN EXPECT + with pytest.raises(sc.UnitError, match='dimensionless argument'): + _ = expr.output_unit + + def test_output_unit_symbolic_power_of_dimensional_base_raises(self): + # WHEN + with pytest.warns(UserWarning, match='not unit-consistent'): + expr = ExpressionComponent( + 'sigma**n * x', + parameters={'sigma': 0.5, 'n': 2.0}, + parameter_units={'sigma': 'meV'}, + ) + + # THEN EXPECT + with pytest.raises(sc.UnitError, match='symbolic power'): + _ = expr.output_unit + + def test_output_unit_exponent_with_unit_raises(self): + # WHEN: the exponent itself carries a unit + with pytest.warns(UserWarning, match='not unit-consistent'): + expr = ExpressionComponent( + 'x**A', parameters={'A': 1.0}, parameter_units={'A': 'meV'}, x_unit='meV' + ) + + # THEN EXPECT + with pytest.raises(sc.UnitError, match='exponent must be dimensionless'): + _ = expr.output_unit + + def test_output_unit_quarter_power_of_dimensional_base_raises(self): + # WHEN: only integer and half-integer powers of quantities with units are supported + with pytest.warns(UserWarning, match='not unit-consistent'): + expr = ExpressionComponent( + 'sigma**0.25 * x', parameters={'sigma': 0.5}, parameter_units={'sigma': 'meV'} + ) + + # THEN EXPECT + with pytest.raises(sc.UnitError, match='raised to the power'): + _ = expr.output_unit + + def test_output_unit_with_x_unit_none(self): + # WHEN: without an x_unit, x is treated as dimensionless in the propagation + with pytest.warns(UserWarning, match='does not match'): + expr = ExpressionComponent( + 'A * x', parameters={'A': 1.0}, parameter_units={'A': 'meV'}, x_unit=None + ) + + # THEN EXPECT + assert sc.Unit(expr.output_unit) == sc.Unit('meV') + + def test_evaluate_raises_for_genuinely_different_x_unit(self): + # WHEN: x carries a different (but compatible) unit than x_unit + expr = ExpressionComponent( + 'A * x', parameters={'A': 1.0}, parameter_units={'A': '1/meV'}, x_unit='meV' + ) + x = sc.linspace('energy', -1.0, 1.0, 11, unit='ueV') + + # THEN EXPECT: ExpressionComponent cannot auto-convert, so evaluating with + # wrongly-scaled x raises instead of silently producing wrong values + with pytest.raises(sc.UnitError, match='cannot auto-convert'): + expr.evaluate(x) + + def test_evaluate_does_not_warn_for_equivalent_unit_spelling(self): + # WHEN: 'ueV' and scipp's canonical micro-sign spelling are the same unit + expr = ExpressionComponent( + 'A * x', parameters={'A': 1.0}, parameter_units={'A': '1/ueV'}, x_unit='ueV' + ) + x = sc.linspace('energy', -1.0, 1.0, 11, unit='\u00b5eV') + + # THEN EXPECT: no spurious warning from the differing spellings + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr.evaluate(x) + + def test_set_unit_warns_when_breaking_consistency(self): + # WHEN: a consistent expression + expr = ExpressionComponent( + 'A * (x - x0)', + parameters={'A': 1.0, 'x0': 0.5}, + parameter_units={'A': '1/meV', 'x0': 'meV'}, + x_unit='meV', + ) + + # THEN EXPECT: relabelling A breaks the output unit + with pytest.warns(UserWarning, match='does not match'): + expr.set_unit('A', '1/ueV') + + +class TestExpressionComponentPhysicalConstants: + def test_kb_constant_value_and_unit(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT: Boltzmann constant in meV/K, a read-only DescriptorNumber + assert expr.kb.value == pytest.approx(0.08617333, rel=1e-6) + assert sc.Unit(str(expr.kb.unit)) == sc.Unit('meV/K') + assert isinstance(expr.kb, DescriptorNumber) + assert not isinstance(expr.kb, Parameter) + + def test_hbar_constant_value_and_unit(self): + # WHEN + expr = ExpressionComponent( + 'exp(-(x * tau / hbar)**2)', parameters={'tau': 1.0}, parameter_units={'tau': 'ps'} + ) + + # THEN EXPECT: hbar in meV*s, a read-only DescriptorNumber + assert expr.hbar.value == pytest.approx(6.582119569509065e-13, rel=1e-6) + assert sc.Unit(str(expr.hbar.unit)) == sc.Unit('meV*s') + assert isinstance(expr.hbar, DescriptorNumber) + assert not isinstance(expr.hbar, Parameter) + + def test_boltzmann_factor_value(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + x = np.array([0.0, 10.0, 25.0]) + + # THEN + result = expr.evaluate(x) + + # EXPECT + expected = np.exp(-x / (0.08617333262145177 * 300.0)) + np.testing.assert_allclose(result, expected, rtol=1e-9) + + def test_constants_are_unit_consistent(self): + # WHEN THEN EXPECT: x/(kb*T) is dimensionless, so no warning is raised + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + assert expr.output_unit == 'dimensionless' + + def test_constant_with_dimensionless_temperature_warns(self): + # WHEN THEN EXPECT: forgetting the unit on T makes x/(kb*T) carry kelvins + with pytest.warns(UserWarning, match='not unit-consistent'): + ExpressionComponent('exp(-x / (kb * T))', parameters={'T': 300.0}) + + def test_constants_property(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT + assert set(expr.constants) == {'kb'} + + def test_constant_not_in_get_all_variables(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT: constants are not fittable variables + names = {p.name for p in expr.get_all_variables()} + assert names == {'T'} + + def test_constant_cannot_be_set(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT + with pytest.raises(AttributeError, match='physical constant'): + expr.kb = 5.0 + + def test_set_unit_on_constant_raises(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT + with pytest.raises(AttributeError, match='physical constant'): + expr.set_unit('kb', 'meV/K') + + def test_parameter_units_for_constant_raises(self): + # WHEN THEN EXPECT + with pytest.raises(ValueError, match='physical constant'): + ExpressionComponent( + 'exp(-x / (kb * T))', + parameters={'T': 300.0}, + parameter_units={'T': 'K', 'kb': 'meV/K'}, + ) + + def test_constant_convert_unit_rescales(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN + expr.kb.convert_unit('eV/K') + + # EXPECT: the value is rescaled with the unit + assert expr.kb.value == pytest.approx(8.617333e-05, rel=1e-6) + + def test_user_parameter_overrides_constant(self): + # WHEN: the user explicitly provides 'kb' as a parameter + with warnings.catch_warnings(): + warnings.simplefilter('error') + expr = ExpressionComponent('exp(-x / (kb * T))', parameters={'kb': 2.0, 'T': 300.0}) + + # THEN EXPECT: it is an ordinary (fittable, dimensionless) parameter + assert not expr.constants + assert expr.kb.value == pytest.approx(2.0) + assert expr.kb.fixed is False + names = {p.name for p in expr.get_all_variables()} + assert names == {'kb', 'T'} + + def test_constant_in_dir(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT + assert 'kb' in dir(expr) + + def test_constant_in_repr(self): + # WHEN + expr = ExpressionComponent( + 'exp(-x / (kb * T))', parameters={'T': 300.0}, parameter_units={'T': 'K'} + ) + + # THEN EXPECT + assert 'kb' in repr(expr)