diff --git a/docs/docs/tutorials/delta_lorentz.ipynb b/docs/docs/tutorials/delta_lorentz.ipynb new file mode 100644 index 00000000..2d0af7b8 --- /dev/null +++ b/docs/docs/tutorials/delta_lorentz.ipynb @@ -0,0 +1,108 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Delta Lorentz\n", + "Model of Delta function and Lorentzian with intensities given by the Debye-Waller factor:\n", + "\n", + "$ I\n", + "= K \\exp \\left( \\frac{-\\langle u^2 \\rangle Q^2}{3} \\right)[A_0 \\delta(E) + (A_1) L(E, \\Gamma)]\n", + "$,\n", + "\n", + "where $K$ is the scale factor, $\\langle u^2 \\rangle$ is the mean square displacement, $Q$ is\n", + "the scattering vector, $A_0$ and $A_1$ are the relative amplitudes of the delta function and\n", + "Lorentzian, respectively, with the constraint that $A_0+A_1=1$, and $L(E, \\Gamma)$ is the\n", + "Lorentzian function with width $\\Gamma$. $A_0$, $A_1$ and the width of the Lorentzian can be\n", + "the same at all $Q$ or be allowed to vary with $Q$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from easydynamics.sample_model.diffusion_model.delta_lorentz import DeltaLorentz\n", + "\n", + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d6e1984", + "metadata": {}, + "outputs": [], + "source": [ + "Q = np.linspace(0.5, 2, 7)\n", + "energy = np.linspace(-2, 2, 501)\n", + "scale = 1.0\n", + "mean_u_squared = 0.5\n", + "A_0 = 0.2\n", + "lorentzian_width = 0.2\n", + "\n", + "diffusion_model = DeltaLorentz(\n", + " scale=scale,\n", + " mean_u_squared=mean_u_squared,\n", + " A_0=A_0,\n", + " lorentzian_width=lorentzian_width,\n", + " allow_Q_variation={'A_0': True},\n", + " Q=Q,\n", + " lorentzian_name='Lorentzian',\n", + " delta_name='Delta function',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "389aa3cd", + "metadata": {}, + "outputs": [], + "source": [ + "component_collections = diffusion_model.get_component_collections()\n", + "cmap = plt.cm.jet\n", + "nQ = len(component_collections)\n", + "plt.figure()\n", + "for Q_index in range(len(component_collections)):\n", + " color = cmap(Q_index / (nQ - 1))\n", + " y = component_collections[Q_index].evaluate(energy)\n", + " plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color)\n", + "\n", + "plt.legend()\n", + "plt.show()\n", + "plt.xlabel('Energy (meV)')\n", + "plt.ylabel('Intensity (arb. units)')\n", + "plt.title('Delta-Lorentz Model')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/docs/tutorials/diffusion_model.ipynb b/docs/docs/tutorials/diffusion_model.ipynb index f3d1571b..148f37c3 100644 --- a/docs/docs/tutorials/diffusion_model.ipynb +++ b/docs/docs/tutorials/diffusion_model.ipynb @@ -6,7 +6,7 @@ "metadata": {}, "source": [ "# Diffusion Model\n", - "We support several standard models of diffusion. Here we show an example of Browniand Translational Diffusion, where the scattering is a Lorentzian with width ($\\Gamma$) given by $\\Gamma = D Q^2$, where $D$ is the diffusion coefficient (in m$^2$/s) and $Q$ is the momentum transfer." + "We support several standard models of diffusion. Here we show an example of Brownian Translational Diffusion, where the scattering is a Lorentzian with width ($\\Gamma$) given by $\\Gamma = D Q^2$, where $D$ is the diffusion coefficient (in m$^2$/s) and $Q$ is the momentum transfer." ] }, { @@ -41,12 +41,10 @@ "diffusion_coefficient = 2.4e-9 # m^2/s\n", "\n", "diffusion_model = BrownianTranslationalDiffusion(\n", - " display_name='DiffusionModel',\n", - " scale=scale,\n", - " diffusion_coefficient=diffusion_coefficient,\n", + " display_name='DiffusionModel', scale=scale, diffusion_coefficient=diffusion_coefficient, Q=Q\n", ")\n", "\n", - "component_collections = diffusion_model.create_component_collections(Q)\n", + "component_collections = diffusion_model.get_component_collections()\n", "\n", "\n", "cmap = plt.cm.jet\n", @@ -87,7 +85,7 @@ ], "metadata": { "kernelspec": { - "display_name": "easydynamics_newbase", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -101,7 +99,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.12" + "version": "3.14.4" } }, "nbformat": 4, diff --git a/docs/docs/tutorials/index.md b/docs/docs/tutorials/index.md index 6c4158a3..0617edb2 100644 --- a/docs/docs/tutorials/index.md +++ b/docs/docs/tutorials/index.md @@ -50,6 +50,9 @@ tutorials. balancing. - [Diffusion model](diffusion_model.ipynb) – Learn how to create and use a model of diffusion. +- [DeltaLorentz](delta_lorentz.ipynb) – Learn how to create and use a + model with a Delta function and a Lorentzian that have a shared + Debye-Waller-like Q-dependence. - [Sample model](sample_model.ipynb) – Learn how to create a model of the scattering from your sample including model components and diffusion models. diff --git a/docs/docs/tutorials/tutorial0_more_advanced.ipynb b/docs/docs/tutorials/tutorial0_more_advanced.ipynb index 2ebc7120..25088780 100644 --- a/docs/docs/tutorials/tutorial0_more_advanced.ipynb +++ b/docs/docs/tutorials/tutorial0_more_advanced.ipynb @@ -33,7 +33,7 @@ "id": "4c8e97b7", "metadata": {}, "source": [ - "As before, we createa an `Experiment` to hold the data, and rebin it." + "As before, we create an `Experiment` to hold the data, and rebin it." ] }, { diff --git a/src/easydynamics/analysis/fit_binding.py b/src/easydynamics/analysis/fit_binding.py index ad11ff60..cbd108dc 100644 --- a/src/easydynamics/analysis/fit_binding.py +++ b/src/easydynamics/analysis/fit_binding.py @@ -257,6 +257,11 @@ def get_parameter_names(self) -> list[str]: modes = self._get_modes() if isinstance(self.model, DiffusionModelBase): + # This needs to be generalised. + # TODO: Generalise this for different diffusion models and modes. # noqa TD002 TD003 + if 'delta' in modes: + return [f'{self.parameter_name} area' for mode in modes] + return [f'{self.parameter_name} {mode}' for mode in modes] return [self.parameter_name] @@ -292,6 +297,9 @@ def _build_diffusion_callable(self, mode: str) -> Callable: if mode == 'width': return lambda x, **_: model.calculate_width(x) + if mode == 'delta': + return lambda x, **_: model.calculate_EISF(x) * model.scale.value + raise ValueError(f'Unknown diffusion mode: {mode}') def _get_modes(self) -> list[str]: diff --git a/src/easydynamics/sample_model/__init__.py b/src/easydynamics/sample_model/__init__.py index 6f75bafd..07f434a1 100644 --- a/src/easydynamics/sample_model/__init__.py +++ b/src/easydynamics/sample_model/__init__.py @@ -16,6 +16,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.sample_model.instrument_model import InstrumentModel from easydynamics.sample_model.resolution_model import ResolutionModel from easydynamics.sample_model.sample_model import SampleModel @@ -26,10 +30,12 @@ 'ComponentCollection', 'DampedHarmonicOscillator', 'DeltaFunction', + 'DeltaLorentz', 'Exponential', 'ExpressionComponent', 'Gaussian', 'InstrumentModel', + 'JumpTranslationalDiffusion', 'Lorentzian', 'Polynomial', 'ResolutionModel', diff --git a/src/easydynamics/sample_model/component_collection.py b/src/easydynamics/sample_model/component_collection.py index 11e6574a..f0bbd2d7 100644 --- a/src/easydynamics/sample_model/component_collection.py +++ b/src/easydynamics/sample_model/component_collection.py @@ -76,7 +76,7 @@ def __init__( Raises ------ TypeError - If unit is not a string or sc.Unit, or if components is not a list of ModelComponent. + If components is not a list of ModelComponent. """ if components is None: components = [] @@ -340,29 +340,6 @@ def free_all_parameters(self) -> None: # Dunder methods # ------------------------------------------------------------------ - def __contains__(self, item: str | ModelComponent) -> bool: - """ - Check if a component with the given name or instance exists in the ComponentCollection. - - Parameters - ---------- - item : str | ModelComponent - The component name or instance to check for. - - Returns - ------- - bool - True if the component exists, False otherwise. - """ - - if isinstance(item, str): - # Check by component name - return any(comp.name == item for comp in self) - if isinstance(item, ModelComponent): - # Check by component instance - return any(comp is item for comp in self) - return False - def __repr__(self) -> str: """ Return a string representation of the ComponentCollection. @@ -374,7 +351,10 @@ def __repr__(self) -> str: """ comp_names = ', '.join(c.name for c in self) or 'No components' - return f"" + return ( + f"ComponentCollection(name='{self.name}', unit='{self.unit}', \n" + f'Components: {comp_names})' + ) def to_dict(self) -> dict: return { 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 0883db57..11f10b1b 100644 --- a/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py +++ b/src/easydynamics/sample_model/diffusion_model/brownian_translational_diffusion.py @@ -12,7 +12,6 @@ from easydynamics.sample_model.diffusion_model.diffusion_model_base import DiffusionModelBase 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 angstrom from easydynamics.utils.utils import hbar @@ -42,9 +41,12 @@ def __init__( self, scale: Numeric = 1.0, diffusion_coefficient: Numeric = 1.0, + Q: Q_type | None = None, unit: str | sc.Unit = 'meV', name: str = 'BrownianTranslationalDiffusion', display_name: str | None = 'BrownianTranslationalDiffusion', + lorentzian_name: str | None = None, + lorentzian_display_name: str | None = None, unique_name: str | None = None, ) -> None: """ @@ -56,12 +58,20 @@ def __init__( Scale factor for the diffusion model. Must be a non-negative number. diffusion_coefficient : Numeric, default=1.0 Diffusion coefficient D in m^2/s. + Q : Q_type | None, default=None + Q values for the model. If None, Q is not set. unit : str | sc.Unit, default='meV' Unit of the diffusion model. Must be convertible to meV. name : str, default='BrownianTranslationalDiffusion' Name of the diffusion model. display_name : str | None, default='BrownianTranslationalDiffusion' 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. + lorentzian_display_name : str | None, default=None + Display name of the Lorentzian component. If None, it will be set to the + lorentzian_name. unique_name : str | None, default=None Unique name of the diffusion model. If None, a unique name will be generated. By default, None. @@ -70,13 +80,27 @@ def __init__( ------ TypeError If scale or diffusion_coefficient is not a number. + + ValueError + If scale or diffusion_coefficient is negative. """ - if not isinstance(scale, Numeric): - raise TypeError('scale must be a number.') + super().__init__( + Q=Q, + unit=unit, + scale=scale, + name=name, + display_name=display_name, + unique_name=unique_name, + lorentzian_name=lorentzian_name, + lorentzian_display_name=lorentzian_display_name, + ) if not isinstance(diffusion_coefficient, Numeric): raise TypeError('diffusion_coefficient must be a number.') + if float(diffusion_coefficient) < 0: + raise ValueError('diffusion_coefficient must be non-negative.') + diffusion_coefficient = Parameter( name='diffusion_coefficient', value=float(diffusion_coefficient), @@ -84,17 +108,13 @@ def __init__( unit='m**2/s', min=0.0, ) - super().__init__( - unit=unit, - scale=scale, - name=name, - display_name=display_name, - unique_name=unique_name, - ) + self._hbar = hbar self._angstrom = angstrom self._diffusion_coefficient = diffusion_coefficient + self._component_collections = self.create_component_collections() + # ------------------------------------------------------------------ # Properties # ------------------------------------------------------------------ @@ -139,13 +159,13 @@ def diffusion_coefficient(self, diffusion_coefficient: Numeric) -> None: # Other methods # ------------------------------------------------------------------ - def calculate_width(self, Q: Q_type) -> np.ndarray: + def calculate_width(self, Q: Q_type | None = None) -> np.ndarray: """ Calculate the half-width at half-maximum (HWHM) for the diffusion model. Parameters ---------- - Q : Q_type + Q : Q_type | None, default=None Scattering vector in 1/angstrom. Returns @@ -153,21 +173,20 @@ def calculate_width(self, Q: Q_type) -> np.ndarray: np.ndarray HWHM values in the unit of the model (e.g., meV). """ - - Q = _validate_and_convert_Q(Q) + Q = self._ensure_Q(Q) unit_conversion_factor = self._hbar * self.diffusion_coefficient / (self._angstrom**2) unit_conversion_factor.convert_unit(self.unit) return Q**2 * unit_conversion_factor.value - def calculate_EISF(self, Q: Q_type) -> np.ndarray: + def calculate_EISF(self, Q: Q_type | None = None) -> np.ndarray: """ Calculate the Elastic Incoherent Structure Factor (EISF) for the Brownian translational diffusion model. Parameters ---------- - Q : Q_type + Q : Q_type | None, default=None Scattering vector in 1/angstrom. Returns @@ -175,16 +194,17 @@ def calculate_EISF(self, Q: Q_type) -> np.ndarray: np.ndarray EISF values (dimensionless). """ - Q = _validate_and_convert_Q(Q) + Q = self._ensure_Q(Q) + return np.zeros_like(Q) - def calculate_QISF(self, Q: Q_type) -> np.ndarray: + def calculate_QISF(self, Q: Q_type | None = None) -> np.ndarray: """ Calculate the Quasi-Elastic Incoherent Structure Factor (QISF). Parameters ---------- - Q : Q_type + Q : Q_type | None, default=None Scattering vector in 1/angstrom. Returns @@ -192,34 +212,17 @@ def calculate_QISF(self, Q: Q_type) -> np.ndarray: np.ndarray QISF values (dimensionless). """ + Q = self._ensure_Q(Q) - Q = _validate_and_convert_Q(Q) return np.ones_like(Q) def create_component_collections( self, - Q: Q_type, - component_name: str = 'Brownian diffusion', - component_display_name: str | None = None, ) -> list[ComponentCollection]: r""" Create ComponentCollection components for the Brownian translational diffusion model at given Q values. - Parameters - ---------- - Q : Q_type - Scattering vector values. - component_name : str, default='Brownian diffusion' - Name of the Brownian diffusion component. - component_display_name : str | None, default=None - Display name of the Brownian diffusion component. - - Raises - ------ - TypeError - If component_display_name is not a string. If component_name is not a string. - Returns ------- list[ComponentCollection] @@ -227,16 +230,10 @@ def create_component_collections( Lorentzian has a width given by $D*Q^2$ and an area given by the scale parameter multiplied by the QISF (which is 1 for this model). """ - Q = _validate_and_convert_Q(Q) - - if not isinstance(component_name, str): - raise TypeError('component_name must be a string.') - - if component_display_name is None: - component_display_name = component_name - - if not isinstance(component_display_name, str): - raise TypeError('component_display_name must be a string.') + Q = self.Q + if Q is None: + self._component_collections = [] + return self._component_collections component_collection_list = [None] * len(Q) # In more complex models, this is used to scale the area of the @@ -254,8 +251,8 @@ def create_component_collections( ) lorentzian_component = Lorentzian( - name=component_name, - display_name=component_display_name, + name=self.lorentzian_name, + display_name=self.lorentzian_display_name, unit=self.unit, ) @@ -284,6 +281,13 @@ def create_component_collections( # Private methods # ------------------------------------------------------------------ + def _on_Q_change(self) -> None: + """ + Update the component collections when Q changes. This is called automatically when the Q + property is set. It regenerates the component collections based on the new Q values. + """ + self._component_collections = self.create_component_collections() + def _write_width_dependency_expression(self, Q: float) -> str: """ Write the dependency expression for the width as a function of Q to make dependent diff --git a/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py b/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py new file mode 100644 index 00000000..d2847d31 --- /dev/null +++ b/src/easydynamics/sample_model/diffusion_model/delta_lorentz.py @@ -0,0 +1,1015 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + + +import numpy as np +import scipp as sc +from easyscience.variable import DescriptorNumber +from easyscience.variable import Parameter + +from easydynamics.sample_model.component_collection import ComponentCollection +from easydynamics.sample_model.components import DeltaFunction +from easydynamics.sample_model.components import Lorentzian +from easydynamics.sample_model.diffusion_model.diffusion_model_base import DiffusionModelBase +from easydynamics.utils.utils import Numeric +from easydynamics.utils.utils import Q_type + +MINIMUM_WIDTH = 1e-10 # To avoid division by zero + + +class DeltaLorentz(DiffusionModelBase): + r""" + Model of Delta function and Lorentzian with intensities given by the Debye-Waller factor. $$ I + = K \exp \left( \frac{-\langle u^2 \rangle Q^2}{3} \right)[A_0 \delta(E) + (A_1) L(E, \Gamma)] + $$, + + where $K$ is the scale factor, $\langle u^2 \rangle$ is the mean square displacement, $Q$ is + the scattering vector, $A_0$ and $A_1$ are the relative amplitudes of the delta function and + Lorentzian, respectively, with the constraint that $A_0+A_1=1$, and $L(E, \Gamma)$ is the + Lorentzian function with width $\Gamma$. $A_0$, $A_1$ and the width of the Lorentzian can be + the same at all $Q$ or be allowed to vary with $Q$. + + + Examples + -------- + >>> Q = np.linspace(0.5, 2, 7) + >>> energy = np.linspace(-2, 2, 501) + >>> scale = 1.0 + >>> mean_u_squared = 0.02 + >>> A_0 = 0.7 + >>> lorentzian_width = 1.0 + >>> model = DeltaLorentz( + ... display_name='DiffusionModel', + ... scale=scale, + ... mean_u_squared=mean_u_squared, + ... A_0=A_0, + ... lorentzian_width=lorentzian_width, + ... allow_Q_variation={'A_0': True, 'lorentzian_width': True}, + ... Q=Q, + ... ) + + See also the tutorials. + """ + + def __init__( + self, + scale: Numeric = 1.0, + mean_u_squared: Numeric = 0.0, + A_0: Numeric = 1.0, + lorentzian_width: Numeric = 1.0, + allow_Q_variation: dict | None = None, + Q: Q_type | None = None, + unit: str | sc.Unit = 'meV', + name: str = 'DeltaLorentz', + display_name: str | None = None, + lorentzian_name: str = 'Lorentzian', + lorentzian_display_name: str | None = None, + delta_name: str = 'Delta function', + delta_display_name: str | None = None, + unique_name: str | None = None, + ) -> None: + """ + Initialize a new DeltaLorentz model. + + Parameters + ---------- + scale : Numeric, default=1.0 + Scale factor for the diffusion model. Must be a non-negative number. + mean_u_squared : Numeric, default=0.0 + Mean square displacement in angstrom^2. + A_0 : Numeric, default=1.0 + Amplitude of the delta function. + lorentzian_width : Numeric, default=1.0 + Width of the Lorentzian function. + allow_Q_variation : dict | None, default=None + Dict describing whether to allow Q variation of A_0 and the Lorentzian width. The dict + can have the keys "A_0" and/or "lorentzian_width", with boolean values indicating + whether to allow Q-dependence for each parameter. If None, no Q-dependence will be + allowed. + Q : Q_type | None, default=None + Q values for the model. If None, Q is not set. + unit : str | sc.Unit, default='meV' + Unit of the diffusion model. Must be convertible to meV. + name : str, default='DeltaLorentz' + Name of the diffusion model. + display_name : str | None, default=None + Display name of the diffusion model. + lorentzian_name : str, default='Lorentzian' + Name of the Lorentzian component. If None, it will be set to the name of the diffusion + model. + 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. + delta_name : str, default='Delta function' + Name of the delta function component. + delta_display_name : str | None, default=None + Display name of the delta function component. If None, it will be set to the display + name of the delta function component. + unique_name : str | None, default=None + Unique name of the diffusion model. If None, a unique name will be generated. By + default, None. + + Raises + ------ + TypeError + If delta_name is not a string or if delta_display_name is not a string or None. + """ + super().__init__( + scale=scale, + unit=unit, + Q=Q, + lorentzian_name=lorentzian_name, + lorentzian_display_name=lorentzian_display_name, + name=name, + display_name=display_name, + unique_name=unique_name, + ) + + # -------------------------------------------------------------- + # Parameters + # -------------------------------------------------------------- + self._mean_u_squared = self._create_mean_u_squared_parameter(mean_u_squared) + + self._A_0, self._A_1 = self._create_A0_A1_parameters(A_0) + + self._lorentzian_width = self._create_lorentzian_width_parameter(lorentzian_width) + + # -------------------------------------------------------------- + # names + # -------------------------------------------------------------- + if not isinstance(delta_name, str): + raise TypeError('delta_name must be a string.') + + if delta_display_name is None: + delta_display_name = delta_name + + if not isinstance(delta_display_name, str): + raise TypeError('delta_display_name must be a string or None.') + + self._delta_name = delta_name + self._delta_display_name = delta_display_name + + # -------------------------------------------------------------- + # Q variation + # -------------------------------------------------------------- + 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() + + # ------------------------------------------------------------------ + # Properties + # ------------------------------------------------------------------ + + @property + def mean_u_squared(self) -> Parameter: + """ + Get the mean square displacement parameter. + + Returns + ------- + Parameter + Mean square displacement in angstrom^2. + """ + return self._mean_u_squared + + @mean_u_squared.setter + def mean_u_squared(self, mean_u_squared: Numeric) -> None: + """ + Set the mean square displacement parameter. + + Parameters + ---------- + mean_u_squared : Numeric + The new value for the mean square displacement in angstrom^2. + + Raises + ------ + TypeError + If mean_u_squared is not a number. + ValueError + If mean_u_squared is negative. + """ + if not isinstance(mean_u_squared, Numeric): + raise TypeError('mean_u_squared must be a number.') + + if float(mean_u_squared) < 0: + raise ValueError('mean_u_squared must be non-negative.') + self._mean_u_squared.value = float(mean_u_squared) + + @property + def A_0(self) -> Parameter: + """ + Get the amplitude of the delta function. + + Returns + ------- + Parameter + Amplitude of the delta function. + """ + return self._A_0 + + @A_0.setter + def A_0(self, A_0: Numeric) -> None: + """ + Set the amplitude of the delta function. + + Parameters + ---------- + A_0 : Numeric + The new value for the amplitude of the delta function. Must be between 0 and 1. + + Raises + ------ + TypeError + If A_0 is not a number. + ValueError + If A_0 is not between 0 and 1. + """ + if not isinstance(A_0, Numeric): + raise TypeError('A_0 must be a number.') + + if not (0 <= float(A_0) <= 1): + raise ValueError('A_0 must be between 0 and 1.') + self._A_0.value = float(A_0) + + @property + def A_1(self) -> Parameter: + """ + Get the amplitude of the Lorentzian function. + + Returns + ------- + Parameter + Amplitude of the Lorentzian function. + """ + return self._A_1 + + @A_1.setter + def A_1(self, _A_1: Numeric) -> None: + """ + A_1 cannot be set directly, as it is a dependent parameter defined as 1 - A_0. To change + A_1, set A_0 to the desired value and A_1 will update accordingly. + + Parameters + ---------- + _A_1 : Numeric + The new value for the amplitude of the Lorentzian function. Is ignored + + Raises + ------ + AttributeError + If an attempt is made to set A_1 directly. + """ + raise AttributeError( + 'A_1 is a dependent parameter and cannot be set directly. ' + 'Set A_0 to change A_1 accordingly.' + ) + + @property + def lorentzian_width(self) -> Parameter: + """ + Get the width of the Lorentzian function. + + Returns + ------- + Parameter + Width of the Lorentzian function. + """ + return self._lorentzian_width + + @lorentzian_width.setter + def lorentzian_width(self, lorentzian_width: Numeric) -> None: + """ + Set the width of the Lorentzian function. + + Parameters + ---------- + lorentzian_width : Numeric + The new value for the width of the Lorentzian function. Must be a non-negative number. + + Raises + ------ + TypeError + If lorentzian_width is not a number. + ValueError + If lorentzian_width is less than the minimum allowed width. + """ + if not isinstance(lorentzian_width, Numeric): + raise TypeError('lorentzian_width must be a number.') + + if float(lorentzian_width) < MINIMUM_WIDTH: + raise ValueError(f'lorentzian_width must be at least {MINIMUM_WIDTH}.') + self._lorentzian_width.value = float(lorentzian_width) + + @property + def delta_name(self) -> str: + """ + Get the name of the delta function component. + + Returns + ------- + str + Name of the delta function component. + """ + return self._delta_name + + @delta_name.setter + def delta_name(self, delta_name: str) -> None: + """ + Set the name of the delta function component. + + Parameters + ---------- + delta_name : str + The new name for the delta function component. + + Raises + ------ + TypeError + If delta_name is not a string. + """ + if not isinstance(delta_name, str): + raise TypeError('delta_name must be a string.') + self._delta_name = delta_name + + @property + def delta_display_name(self) -> str: + """ + Get the display name of the delta function component. + + Returns + ------- + str + Display name of the delta function component. + """ + return self._delta_display_name + + @delta_display_name.setter + def delta_display_name(self, delta_display_name: str | None) -> None: + """ + Set the display name of the delta function component. + + Parameters + ---------- + delta_display_name : str | None + The new display name for the delta function component. + + Raises + ------ + TypeError + If delta_display_name is not a string or None. + """ + if not isinstance(delta_display_name, (str, type(None))): + raise TypeError('delta_display_name must be a string or None.') + self._delta_display_name = delta_display_name + + # ------------------------------------------------------------------ + # Other methods + # ------------------------------------------------------------------ + + def calculate_width(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the half-width at half-maximum (HWHM) for the diffusion model. If the width is + allowed to vary with Q then the Q stored in the model is used and the input is ignored. If + the width is not allowed to vary then the same width is returned for all Q values. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. + + Returns + ------- + np.ndarray + HWHM values in the unit of the model (e.g., meV). + """ + if self._allow_Q_variation['lorentzian_width'] is True: + widths = [lorentzian_width.value for lorentzian_width in self._lorentzian_width_list] + return np.array(widths) + + Q = self._ensure_Q(Q) + + widths = self.lorentzian_width.value * np.ones_like(Q) + + return np.array(widths) + + def calculate_EISF(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Elastic Incoherent Structure Factor (EISF) for the diffusion model. + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. + + Returns + ------- + np.ndarray + EISF values (dimensionless). + """ + Q = self._ensure_Q(Q) + if self._allow_Q_variation['A_0'] is True: + A_0_values = [A_0_.value for A_0_ in self._A_0_list] + return np.exp(-self.mean_u_squared.value * Q**2 / 3) * np.array(A_0_values) + + A_0_values = [self.A_0.value] * len(Q) + return np.exp(-self.mean_u_squared.value * Q**2 / 3) * np.array(A_0_values) + + def calculate_QISF(self, Q: Q_type = None) -> np.ndarray: + """ + Calculate the Quasi-Elastic Incoherent Structure Factor (QISF). + + Parameters + ---------- + Q : Q_type, default=None + Scattering vector in 1/angstrom. + + Returns + ------- + np.ndarray + QISF values (dimensionless). + """ + Q = self._ensure_Q(Q) + if self._allow_Q_variation['A_0'] is True: + A_1_values = [A_1_.value for A_1_ in self._A_1_list] + return np.exp(-self.mean_u_squared.value * Q**2 / 3) * np.array(A_1_values) + + A_1_values = [self.A_1.value] * len(Q) + return np.exp(-self.mean_u_squared.value * Q**2 / 3) * np.array(A_1_values) + + def create_component_collections( + self, + ) -> list[ComponentCollection]: + r""" + Create ComponentCollection components for the DeltaLorentz model at given Q values. + + Returns + ------- + list[ComponentCollection] + List of ComponentCollections with Lorentzian and delta function components for each Q + value. + """ + Q = self.Q + if Q is None: + return [] + + 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 + + if self._allow_Q_variation['lorentzian_width'] is True: + lorentzian_width_list = self._create_lorentzian_width_parameter_list( + self.lorentzian_width + ) + self._lorentzian_width_list = lorentzian_width_list + + component_collection_list = [None] * len(Q) + for i, Q_value in enumerate(Q): + component_collection_list[i] = ComponentCollection( + display_name=f'{self.display_name}_Q{Q_value:.2f}', + unit=self.unit, + ) + + # ------------------------------# + # Create Lorentzian + # ------------------------------# + if self._allow_Q_variation['lorentzian_width'] is True: + width = self._lorentzian_width_list[i] + else: + width = 1.0 + lorentzian_component = Lorentzian( + name=self.lorentzian_name, + display_name=self.lorentzian_display_name, + unit=self.unit, + width=width, + ) + + # If the width is allowed to vary with Q it is independent. + # If the width is not allowed to vary with Q it must be made + # dependent on the width parameter of the model. + if self._allow_Q_variation['lorentzian_width'] is False: + dependency_map = self._write_width_dependency_map_expression() + + lorentzian_component.width.make_dependent_on( + dependency_expression=self._write_lorz_width_dependency_expression(Q_value), + dependency_map=dependency_map, + desired_unit=self.unit, + ) + # The area is always a dependent parameter in this model, as + # it depends on the scale, mean_u_squared and A_1 parameters + # of the model. If A_1 is allowed to vary with Q, the area + # will also depend on the specific A_1 parameter for that Q + # value. If A_1 is not allowed to vary with Q, the area will + # depend on the single A_1 parameter of the model. + if self._allow_Q_variation['A_0'] is True: + dependency_map = self._write_lorz_area_dependency_map_expression(i) + else: + dependency_map = self._write_lorz_area_dependency_map_expression(None) + + lorentzian_component.area.make_dependent_on( + dependency_expression=self._write_lorz_area_dependency_expression(Q_value), + dependency_map=dependency_map, + ) + + component_collection_list[i].append_component(lorentzian_component) + + # ------------------------------# + # Create delta function + # ------------------------------# + + delta_component = DeltaFunction( + name=self.delta_name, + display_name=self.delta_display_name, + unit=self.unit, + ) + + if self._allow_Q_variation['A_0'] is True: + dependency_map = self._write_delta_area_dependency_map_expression(i) + else: + dependency_map = self._write_delta_area_dependency_map_expression(None) + + delta_component.area.make_dependent_on( + dependency_expression=self._write_delta_area_dependency_expression(Q_value), + dependency_map=dependency_map, + ) + + component_collection_list[i].append_component(delta_component) + + return component_collection_list + + def get_all_variables(self, Q_index: int | None = None) -> list[DescriptorNumber]: + """ + Get a list of all variables (Parameters and Descriptors) in the model. + + Parameters + ---------- + Q_index : int | None, default=None + The index of the Q value for which to get the variables. If None, variables for all Q + values will be included. + + Returns + ------- + list[DescriptorNumber] + List of all variables in the model. + """ + + variables = [self.scale, self.mean_u_squared] + if self._allow_Q_variation['A_0'] is True: + if Q_index is None: + variables.extend(self._A_0_list) + variables.extend(self._A_1_list) + else: + variables.append(self._A_0_list[Q_index]) + variables.append(self._A_1_list[Q_index]) + else: + variables.append(self.A_0) + variables.append(self.A_1) + + if self._allow_Q_variation['lorentzian_width'] is True: + if Q_index is None: + variables.extend(self._lorentzian_width_list) + else: + variables.append(self._lorentzian_width_list[Q_index]) + else: + variables.append(self.lorentzian_width) + + return variables + + # ------------------------------------------------------------------ + # Private methods for init + # ------------------------------------------------------------------ + + def _create_Q_variation_dict(self, allow_Q_variation: dict | None) -> dict: + """ + Create a dict for the allow_Q_variation attribute, ensuring that it has the correct keys + and default values. + + Parameters + ---------- + allow_Q_variation : dict | None + Dict describing whether to allow Q variation of A_0 and the Lorentzian width. + + Raises + ------ + TypeError + If allow_Q_variation is not a dict or None. + ValueError + If allow_Q_variation contains unknown keys. + + Returns + ------- + dict + A dict with keys 'A_0' and 'lorentzian_width', indicating whether to allow Q variation + for each parameter. + """ + + allow_Q_variation_default = { + 'A_0': False, + 'lorentzian_width': False, + } + allowed_keys = set(allow_Q_variation_default) + + if allow_Q_variation is None: + allow_Q_variation = {} + if not isinstance(allow_Q_variation, dict): + raise TypeError('allow_Q_variation must be a dict or None.') + + unknown_keys = set(allow_Q_variation) - allowed_keys + if unknown_keys: + raise ValueError(f'Unknown keys in allow_Q_variation: {unknown_keys}') + + return {**allow_Q_variation_default, **allow_Q_variation} + + def _create_mean_u_squared_parameter(self, mean_u_squared: Numeric) -> Parameter: + """ + Create the mean square displacement parameter. + + Parameters + ---------- + mean_u_squared : Numeric + The value for the mean square displacement in angstrom^2. + + Raises + ------ + TypeError + If mean_u_squared is not a number. + ValueError + If mean_u_squared is negative. + + Returns + ------- + Parameter + The created mean square displacement parameter. + """ + + if not isinstance(mean_u_squared, Numeric): + raise TypeError('mean_u_squared must be a number.') + + if float(mean_u_squared) < 0: + raise ValueError('mean_u_squared must be non-negative.') + + return Parameter( + name='mean_u_squared', + value=float(mean_u_squared), + fixed=False, + min=0.0, + unit='angstrom**2', + ) + + def _create_A0_A1_parameters(self, A_0: Numeric) -> tuple[Parameter, Parameter]: + """ + Create the A_0 and A_1 parameters. + + Parameters + ---------- + A_0 : Numeric + The value for the A_0 parameter. + + Raises + ------ + TypeError + If A_0 is not a number. + ValueError + If A_0 is not between 0 and 1. + + Returns + ------- + tuple[Parameter, Parameter] + A tuple containing the A_0 and A_1 parameters. + """ + if not isinstance(A_0, Numeric): + raise TypeError('A_0 must be a number.') + + if float(A_0) < 0 or float(A_0) > 1: + raise ValueError('A_0 must be between 0 and 1.') + + A_0 = Parameter( + name='A_0', + value=float(A_0), + fixed=False, + min=0.0, + max=1.0, + ) + + A_1 = Parameter.from_dependency( + name='A_1', + dependency_expression='1 - A_0', + dependency_map={'A_0': A_0}, + ) + return A_0, A_1 + + def _create_lorentzian_width_parameter(self, lorentzian_width: Numeric) -> Parameter: + """ + Create the Lorentzian width parameter. + + Parameters + ---------- + lorentzian_width : Numeric + The value for the Lorentzian width parameter. + + Raises + ------ + TypeError + If lorentzian_width is not a number. + ValueError + If lorentzian_width is less than the minimum width. + + Returns + ------- + Parameter + The created Lorentzian width parameter. + """ + + if not isinstance(lorentzian_width, Numeric): + raise TypeError('lorentzian_width must be a number.') + + if float(lorentzian_width) < MINIMUM_WIDTH: + raise ValueError(f'lorentzian_width must be at least {MINIMUM_WIDTH}.') + + return Parameter( + name='lorentzian_width', + value=float(lorentzian_width), + fixed=False, + min=MINIMUM_WIDTH, + unit=self.unit, + ) + + def _create_A0_A1_parameter_lists( + self, + A_0: Parameter, + ) -> tuple[list[Parameter], list[Parameter]]: + """ + Create lists of A_0 and A_1 parameters for each Q value. + + Parameters + ---------- + A_0 : Parameter + The A_0 parameter to use as the base for creating the A_0 parameters for each Q value. + + Returns + ------- + tuple[list[Parameter], list[Parameter]] + A tuple containing two lists: the first list contains the A_0 parameters for each Q + value, and the second list contains the A_1 parameters for each Q value. + """ + A_0_list = [] + A_1_list = [] + for _ in self.Q: + a0 = Parameter( + name='A_0', + value=float(A_0.value), + fixed=False, + min=0.0, + max=1.0, + ) + + a1 = Parameter.from_dependency( + name='A_1', + dependency_expression='1 - A_0', + dependency_map={'A_0': a0}, + ) + + A_0_list.append(a0) + A_1_list.append(a1) + + return A_0_list, A_1_list + + def _create_lorentzian_width_parameter_list( + self, + lorentzian_width: Parameter, + ) -> list[Parameter]: + """ + Create a list of Lorentzian width parameters for each Q value. + + Parameters + ---------- + lorentzian_width : Parameter + The Lorentzian width parameter to use as the base for creating the Lorentzian width + parameters for each Q value. + + Returns + ------- + list[Parameter] + A list containing the Lorentzian width parameters for each Q value. + """ + return [ + Parameter( + name=f'{self.lorentzian_name} width', + value=float(lorentzian_width.value), + fixed=False, + min=MINIMUM_WIDTH, + unit=self.unit, + ) + for _ in self.Q + ] + + # ------------------------------------------------------------------ + # Private methods + # ------------------------------------------------------------------ + + 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 = [] + + 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() + + def _write_lorz_width_dependency_expression(self, Q: float) -> str: + """ + Write the dependency expression for the width as a function of Q to make dependent + Parameters. + + Parameters + ---------- + Q : float + Scattering vector in 1/angstrom. + + Raises + ------ + TypeError + If Q is not a float. + + Returns + ------- + str + Dependency expression for the width. + """ + if not isinstance(Q, (float)): + raise TypeError('Q must be a float.') + + return 'lorentzian_width' + + def _write_width_dependency_map_expression( + self, + ) -> dict[str, DescriptorNumber]: + """ + Write the dependency map expression to make dependent Parameters. + + Returns + ------- + dict[str, DescriptorNumber] + Dependency map for the width. + """ + return { + 'lorentzian_width': self.lorentzian_width, + } + + def _write_lorz_area_dependency_expression(self, Q: float) -> str: + """ + Write the dependency expression for the area to make dependent Parameters. + + Parameters + ---------- + Q : float + Scattering vector in 1/angstrom . + + Raises + ------ + TypeError + If Q is not a float. + + Returns + ------- + str + Dependency expression for the area. + """ + if not isinstance(Q, (float)): + raise TypeError('Q must be a float.') + + return f'scale * exp(-mean_u_squared.value * {Q}**2 / 3) * A_1' + + def _write_lorz_area_dependency_map_expression( + self, Q_index: int | None + ) -> dict[str, DescriptorNumber]: + """ + Write the dependency map expression to make dependent Parameters. + + Parameters + ---------- + Q_index : int | None + Index of the Q value for which to write the dependency map. If None, write the + dependency map for the case where A_1 is not Q-dependent. + + Returns + ------- + dict[str, DescriptorNumber] + Dependency map for the area. + """ + if Q_index is None: + return { + 'scale': self.scale, + 'mean_u_squared': self.mean_u_squared, + 'A_1': self.A_1, + } + + return { + 'scale': self.scale, + 'mean_u_squared': self.mean_u_squared, + 'A_1': self._A_1_list[Q_index], + } + + def _write_delta_area_dependency_expression(self, Q: float) -> str: + """ + Write the dependency expression for the area to make dependent Parameters. + + Parameters + ---------- + Q : float + Scattering vector in 1/angstrom. + + Raises + ------ + TypeError + If Q is not a float. + + Returns + ------- + str + Dependency expression for the area. + """ + if not isinstance(Q, (float)): + raise TypeError('Q must be a float.') + + return f'scale * exp(-mean_u_squared.value * {Q}**2 / 3) * A_0' + + def _write_delta_area_dependency_map_expression( + self, + Q_index: int | None, + ) -> dict[str, DescriptorNumber]: + """ + Write the dependency map expression to make dependent Parameters. + + Parameters + ---------- + Q_index : int | None + Index of the Q value for which to write the dependency map. If None, write the + dependency map for the case where A_0 is not Q-dependent. + + Returns + ------- + dict[str, DescriptorNumber] + Dependency map for the area. + """ + if Q_index is None: + return { + 'scale': self.scale, + 'mean_u_squared': self.mean_u_squared, + 'A_0': self.A_0, + } + return { + 'scale': self.scale, + 'mean_u_squared': self.mean_u_squared, + 'A_0': self._A_0_list[Q_index], + } + + # ------------------------------------------------------------------ + # dunder methods + # ------------------------------------------------------------------ + + def __repr__(self) -> str: + """ + String representation of the DeltaLorentz model. + + Returns + ------- + str + String representation of the DeltaLorentz model. + """ + return ( + f'DeltaLorentz(display_name={self.display_name},' + f'unit={self.unit}, \n' + f' mean_u_squared={self.mean_u_squared}, \n' + f' A_0={self.A_0}, A_1={self.A_1}, \n' + f' lorentzian_width={self.lorentzian_width}, \n' + f' scale={self.scale})' + ) 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 95732680..12565986 100644 --- a/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py +++ b/src/easydynamics/sample_model/diffusion_model/diffusion_model_base.py @@ -1,13 +1,17 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import numpy as np import scipp as sc from easyscience.variable import DescriptorNumber from easyscience.variable import Parameter from scipp import UnitError from easydynamics.base_classes.easydynamics_modelbase import EasyDynamicsModelBase +from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.utils.utils import Numeric +from easydynamics.utils.utils import Q_type +from easydynamics.utils.utils import _validate_and_convert_Q class DiffusionModelBase(EasyDynamicsModelBase): @@ -16,9 +20,12 @@ class DiffusionModelBase(EasyDynamicsModelBase): def __init__( self, scale: Numeric = 1.0, + Q: Q_type | None = None, unit: str | sc.Unit = 'meV', name: str = 'DiffusionModel', - display_name: str | None = 'MyDiffusionModel', + display_name: str | None = 'DiffusionModel', + lorentzian_name: str | None = None, + lorentzian_display_name: str | None = None, unique_name: str | None = None, ) -> None: """ @@ -28,12 +35,20 @@ def __init__( ---------- scale : Numeric, default=1.0 Scale factor for the diffusion model. Must be a non-negative number. + Q : Q_type | None, default=None + Q values for the model. If None, Q is not set. unit : str | sc.Unit, default='meV' Unit of the diffusion model. Must be convertible to meV. name : str, default='DiffusionModel' Name of the diffusion model. - display_name : str | None, default='MyDiffusionModel' + display_name : str | None, default='DiffusionModel' 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. + lorentzian_display_name : str | None, default=None + Display name of the Lorentzian component. If None, it will be set to the + lorentzian_name. unique_name : str | None, default=None Unique name of the diffusion model. If None, a unique name will be generated. By default, None. @@ -44,8 +59,12 @@ def __init__( If scale is not a number. UnitError If unit is not a string or scipp Unit, or if it cannot be converted to meV. + ValueError + If scale is negative. """ + self._Q = _validate_and_convert_Q(Q) + try: test = DescriptorNumber(name='test', value=1, unit=unit) test.convert_unit('meV') @@ -57,10 +76,33 @@ def __init__( if not isinstance(scale, Numeric): raise TypeError('scale must be a number.') + if float(scale) < 0: + raise ValueError('scale must be non-negative.') + scale = Parameter(name='scale', value=float(scale), fixed=False, min=0.0, unit=unit) + self._scale = scale super().__init__(unit=unit, name=name, display_name=display_name, unique_name=unique_name) - self._scale = scale + + if lorentzian_name is None: + lorentzian_name = name + + if not isinstance(lorentzian_name, str): + raise TypeError('lorentzian_name must be a string.') + + if lorentzian_display_name is None: + lorentzian_display_name = lorentzian_name + + if not isinstance(lorentzian_display_name, str): + raise TypeError('lorentzian_display_name must be a string or None.') + + self._lorentzian_name = lorentzian_name + self._lorentzian_display_name = lorentzian_display_name + + if self.Q is None: + self._component_collections = [] + else: + self._component_collections = [ComponentCollection()] * len(self.Q) # ------------------------------------------------------------------ # Properties @@ -102,6 +144,248 @@ def scale(self, scale: Numeric) -> None: raise ValueError('scale must be non-negative.') self._scale.value = float(scale) + @property + def Q(self) -> np.ndarray | None: + """ + Get the Q values of the SampleModel. + + Returns + ------- + np.ndarray | None + The Q values of the SampleModel, or None if not set. + """ + return self._Q + + @Q.setter + def Q(self, value: Q_type | None) -> None: + """ + Set the Q values of the SampleModel. + + If Q is already set, it throws an error if the new Q values are not similar to the old + ones. To change Q values, first run clear_Q(). + + Parameters + ---------- + value : Q_type | None + The new Q values to set. If None, Q values are not changed. + + Raises + ------ + ValueError + If the new Q values are not similar to the old ones when Q is already set. + """ + if value is None: + return + old_Q = self._Q + new_Q = _validate_and_convert_Q(value) + + if old_Q is None: + self._Q = new_Q + self._on_Q_change() + return + + if len(old_Q) != len(new_Q) or not np.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().' + ) + + @property + def lorentzian_name(self) -> str: + """ + Get the name of the Lorentzian component. + + Returns + ------- + str + Name of the Lorentzian component. + """ + return self._lorentzian_name + + @lorentzian_name.setter + def lorentzian_name(self, lorentzian_name: str) -> None: + """ + Set the name of the Lorentzian component. + + Parameters + ---------- + lorentzian_name : str + The new name for the Lorentzian component. + + Raises + ------ + TypeError + If lorentzian_name is not a string. + """ + if not isinstance(lorentzian_name, str): + raise TypeError('lorentzian_name must be a string.') + self._lorentzian_name = lorentzian_name + + @property + def lorentzian_display_name(self) -> str | None: + """ + Get the display name of the Lorentzian component. + + Returns + ------- + str | None + Display name of the Lorentzian component, or None if not set. + """ + return self._lorentzian_display_name + + @lorentzian_display_name.setter + def lorentzian_display_name(self, lorentzian_display_name: str | None) -> None: + """ + Set the display name of the Lorentzian component. + + Parameters + ---------- + lorentzian_display_name : str | None + The new display name for the Lorentzian component. + + Raises + ------ + TypeError + If lorentzian_display_name is not a string or None. + """ + if not isinstance(lorentzian_display_name, (str, type(None))): + raise TypeError('lorentzian_display_name must be a string or None.') + self._lorentzian_display_name = lorentzian_display_name + + def clear_Q(self, confirm: bool = False) -> None: + """ + Clear the Q values of the SampleModel, removing all component collections and their + associated Parameters. + + Parameters + ---------- + confirm : bool, default=False + Confirmation to clear Q values. + + Raises + ------ + ValueError + If confirm is not True. + """ + if not confirm: + raise ValueError( + 'Clearing Q values requires confirmation. Set confirm=True to proceed.' + ) + self._Q = None + self._on_Q_change() + + # ------------------------------------------------------------------ + # private methods + # ------------------------------------------------------------------ + + def create_component_collections(self) -> list[ComponentCollection]: + """ + Create the ComponentCollections for the diffusion model based on the current Q values. + + Returns + ------- + list[ComponentCollection] + A list of ComponentCollections corresponding to the current Q values. + """ + if self.Q is None: + self._component_collections = [] + return self._component_collections + + self._component_collections = [ComponentCollection()] * len(self.Q) + + return self._component_collections + + def get_component_collections( + self, Q_index: int | None = None + ) -> ComponentCollection | list[ComponentCollection]: + """ + Get the ComponentCollection at the given Q index. + + Parameters + ---------- + Q_index : int | None, default=None + The index of the desired ComponentCollection. If None, all ComponentCollections are + returned. + + Raises + ------ + TypeError + If Q_index is not an int. + IndexError + If Q_index is out of bounds for the number of ComponentCollections. + + Returns + ------- + ComponentCollection | list[ComponentCollection] + The ComponentCollection at the specified Q index. If Q_index is None, a list of all + ComponentCollections is returned. + """ + if Q_index is None: + return self._component_collections + + if not isinstance(Q_index, int): + raise TypeError(f'Q_index must be an int, got {type(Q_index).__name__}') + if Q_index < 0 or Q_index >= len(self._component_collections): + raise IndexError( + f'Q_index {Q_index} is out of bounds for component collections ' + f'of length {len(self._component_collections)}' + ) + return self._component_collections[Q_index] + + def get_all_variables( + self, + Q_index: int | None = None, # noqa ARG002 + ) -> list[Parameter]: + """ + Get all Parameters from the ComponentCollection at the given Q index. + + Parameters + ---------- + Q_index : int | None, default=None + Unused parameter for compatibility with SampleModel. + + + Returns + ------- + list[Parameter] + A list of all Parameters from the specified ComponentCollection(s). + """ + return super().get_all_variables() + + # ------------------------------------------------------------------ + # private methods + # ------------------------------------------------------------------ + + def _on_Q_change(self) -> None: + """Handle changes to the Q values.""" + + def _ensure_Q(self, Q: Q_type) -> np.ndarray: + """ + Convert Q to a numpy array, ensuring it is not None. Uses the stored Q if no input is + given. + + Parameters + ---------- + Q : Q_type + The Q to be checked + + Returns + ------- + np.ndarray + The validated and converted Q values. + + Raises + ------ + ValueError + If the provided Q and self.Q are both None + """ + if Q is None: + Q = self.Q + if Q is None: + raise ValueError('Q must be provided either as an argument or set in the model.') + + return _validate_and_convert_Q(Q) + # ------------------------------------------------------------------ # dunder methods # ------------------------------------------------------------------ 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 0d4b4441..160f1476 100644 --- a/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py +++ b/src/easydynamics/sample_model/diffusion_model/jump_translational_diffusion.py @@ -12,7 +12,6 @@ from easydynamics.sample_model.diffusion_model.diffusion_model_base import DiffusionModelBase 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 angstrom from easydynamics.utils.utils import hbar @@ -51,9 +50,12 @@ def __init__( scale: Numeric = 1.0, diffusion_coefficient: Numeric = 1.0, relaxation_time: Numeric = 1.0, + Q: Q_type | None = None, unit: str | sc.Unit = 'meV', name: str = 'JumpTranslationalDiffusion', display_name: str | None = 'JumpTranslationalDiffusion', + lorentzian_name: str | None = None, + lorentzian_display_name: str | None = None, unique_name: str | None = None, ) -> None: """ @@ -67,12 +69,20 @@ def __init__( Diffusion coefficient D in m^2/s. relaxation_time : Numeric, default=1.0 Relaxation time t in ps. + Q : Q_type | None, default=None + Q values for the model. If None, Q is not set. unit : str | sc.Unit, default='meV' Unit of the diffusion model. Must be convertible to meV. name : str, default='JumpTranslationalDiffusion' Name of the diffusion model. display_name : str | None, default='JumpTranslationalDiffusion' 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. + 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 unique_name : str | None, default=None Unique name of the diffusion model. If None, a unique name will be generated. By default, None. @@ -81,26 +91,38 @@ def __init__( ------ TypeError If scale, diffusion_coefficient, or relaxation_time are not numbers. + ValueError + If scale, diffusion_coefficient, or relaxation_time are negative. """ super().__init__( + Q=Q, + unit=unit, + scale=scale, name=name, display_name=display_name, + lorentzian_name=lorentzian_name, + lorentzian_display_name=lorentzian_display_name, unique_name=unique_name, - unit=unit, - scale=scale, ) if not isinstance(diffusion_coefficient, Numeric): raise TypeError('diffusion_coefficient must be a number.') + if float(diffusion_coefficient) < 0: + raise ValueError('diffusion_coefficient must be non-negative.') + if not isinstance(relaxation_time, Numeric): raise TypeError('relaxation_time must be a number.') + if float(relaxation_time) < 0: + raise ValueError('relaxation_time must be non-negative.') + diffusion_coefficient = Parameter( name='diffusion_coefficient', value=float(diffusion_coefficient), fixed=False, unit='m**2/s', + min=0.0, ) relaxation_time = Parameter( @@ -108,6 +130,7 @@ def __init__( value=float(relaxation_time), fixed=False, unit='ps', + min=0.0, ) self._hbar = hbar @@ -115,6 +138,8 @@ def __init__( self._diffusion_coefficient = diffusion_coefficient self._relaxation_time = relaxation_time + self._component_collections = self.create_component_collections() + ################################ # Properties ################################ @@ -194,15 +219,16 @@ def relaxation_time(self, relaxation_time: Numeric) -> None: # Other methods ################################ - def calculate_width(self, Q: Q_type) -> np.ndarray: + def calculate_width(self, Q: Q_type | None = None) -> np.ndarray: r""" Calculate the half-width at half-maximum (HWHM) for the diffusion model. $\Gamma(Q) = Q^2/(1+D t Q^2)$, where $D$ is the diffusion coefficient and $t$ is the relaxation time. Parameters ---------- - Q : Q_type - Scattering vector in 1/angstrom. Can be a single value or an array of values. + Q : Q_type | None, default=None + Scattering vector in 1/angstrom. Can be a single value or an array of values. If None, + Q values stored in the model are used. Returns ------- @@ -210,7 +236,7 @@ def calculate_width(self, Q: Q_type) -> np.ndarray: HWHM values in the unit of the model (e.g., meV). """ - Q = _validate_and_convert_Q(Q) + Q = self._ensure_Q(Q) unit_conversion_factor_numerator = ( self._hbar * self.diffusion_coefficient / (self._angstrom**2) @@ -242,7 +268,8 @@ def calculate_EISF(self, Q: Q_type) -> np.ndarray: np.ndarray EISF values (dimensionless). """ - Q = _validate_and_convert_Q(Q) + Q = self._ensure_Q(Q) + return np.zeros_like(Q) def calculate_QISF(self, Q: Q_type) -> np.ndarray: @@ -259,44 +286,27 @@ def calculate_QISF(self, Q: Q_type) -> np.ndarray: np.ndarray QISF values (dimensionless). """ - Q = _validate_and_convert_Q(Q) + Q = self._ensure_Q(Q) + return np.ones_like(Q) def create_component_collections( self, - Q: Q_type, - component_name: str = 'Jump translational diffusion', - component_display_name: str = 'Jump translational diffusion', ) -> list[ComponentCollection]: """ Create ComponentCollection components for the diffusion model at given Q values. - Parameters - ---------- - Q : Q_type - Scattering vector in 1/angstrom. Can be a single value or an array of values. - component_name : str, default='Jump translational diffusion' - Name of the Jump Diffusion Lorentzian component. - component_display_name : str, default='Jump translational diffusion' - Name of the Jump Diffusion Lorentzian component. - - Raises - ------ - TypeError - If component_display_name is not a string. If component_name is not a string. + TypeError If component_display_name is not a string. If component_name is not a string. Returns ------- list[ComponentCollection] List of ComponentCollections with Jump Diffusion Lorentzian components. """ - Q = _validate_and_convert_Q(Q) - - if not isinstance(component_display_name, str): - raise TypeError('component_display_name must be a string.') - - if not isinstance(component_name, str): - raise TypeError('component_name must be a string.') + Q = self.Q + if Q is None: + self._component_collections = [] + return self._component_collections component_collection_list = [None] * len(Q) # In more complex models, this is used to scale the area of the @@ -314,8 +324,8 @@ def create_component_collections( ) lorentzian_component = Lorentzian( - name=component_name, - display_name=component_display_name, + name=self.lorentzian_name, + display_name=self.lorentzian_display_name, unit=self.unit, ) @@ -343,6 +353,12 @@ def create_component_collections( ################################ # Private methods ################################ + def _on_Q_change(self) -> None: + """ + Update the component collections when Q changes. This is called automatically when the Q + property is set. It regenerates the component collections based on the new Q values. + """ + self._component_collections = self.create_component_collections() def _write_width_dependency_expression(self, Q: float) -> str: """ @@ -380,8 +396,8 @@ def _write_width_dependency_map_expression(self) -> dict[str, DescriptorNumber]: Dependency map for the width. """ return { - 'D': self._diffusion_coefficient, - 't': self._relaxation_time, + 'D': self.diffusion_coefficient, + 't': self.relaxation_time, 'hbar': self._hbar, 'angstrom': self._angstrom, } @@ -421,7 +437,7 @@ def _write_area_dependency_map_expression(self) -> dict[str, DescriptorNumber]: Dependency map for the area. """ return { - 'scale': self._scale, + 'scale': self.scale, } ################################ @@ -439,6 +455,6 @@ def __repr__(self) -> str: """ return ( f'JumpTranslationalDiffusion(name={self.name}, display_name={self.display_name},\n ' - f' diffusion_coefficient={self._diffusion_coefficient}, \n' - f' scale={self._scale})' + f' diffusion_coefficient={self.diffusion_coefficient}, \n' + f' scale={self.scale})' ) diff --git a/src/easydynamics/sample_model/model_base.py b/src/easydynamics/sample_model/model_base.py index 4895dc6b..aa76b384 100644 --- a/src/easydynamics/sample_model/model_base.py +++ b/src/easydynamics/sample_model/model_base.py @@ -141,39 +141,6 @@ def clear_components(self) -> None: # Properties # ------------------------------------------------------------------ - @property - def unit(self) -> str | sc.Unit | None: - """ - Get the unit of the SampleModel. - - Returns - ------- - str | sc.Unit | None - The unit of the SampleModel. - """ - - return self._unit - - @unit.setter - def unit(self, _unit_str: str) -> None: - """ - Unit is read-only and cannot be set directly. - - Parameters - ---------- - _unit_str : str - The new unit to set (ignored). - - Raises - ------ - AttributeError - Always raised to indicate that the unit is read-only. - """ - raise AttributeError( - f'Unit is read-only. Use convert_unit to change the unit between allowed types ' - f'or create a new {self.__class__.__name__} with the desired unit.' - ) - @property def components(self) -> list[ModelComponent]: """ @@ -410,12 +377,12 @@ def normalize_area(self) -> None: def _generate_component_collections(self) -> None: """Generate ComponentCollections for each Q value.""" - if self._Q is None: + if self.Q is None: self._component_collections = [] return self._component_collections = [] - for _ in self._Q: + for _ in self.Q: self._component_collections.append(copy(self._components)) def _on_Q_change(self) -> None: diff --git a/src/easydynamics/sample_model/sample_model.py b/src/easydynamics/sample_model/sample_model.py index bcb598f7..83325506 100644 --- a/src/easydynamics/sample_model/sample_model.py +++ b/src/easydynamics/sample_model/sample_model.py @@ -15,6 +15,7 @@ from easydynamics.utils import detailed_balance_factor from easydynamics.utils.utils import Numeric from easydynamics.utils.utils import Q_type +from easydynamics.utils.utils import _validate_and_convert_Q class SampleModel(ModelBase): @@ -73,6 +74,7 @@ def __init__( ValueError If temperature is negative. """ + if diffusion_models is None: self._diffusion_models = [] elif isinstance(diffusion_models, DiffusionModelBase): @@ -87,6 +89,10 @@ 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, @@ -142,7 +148,7 @@ def append_diffusion_model(self, diffusion_model: DiffusionModelBase) -> None: raise TypeError( f'diffusion_model must be a DiffusionModelBase, got {type(diffusion_model).__name__}' # noqa: E501 ) - + diffusion_model.Q = self.Q self._diffusion_models.append(diffusion_model) self._generate_component_collections() @@ -160,19 +166,19 @@ def remove_diffusion_model(self, name: str) -> None: ValueError If no DiffusionModel with the given name is found. """ - for i, dm in enumerate(self._diffusion_models): + for i, dm in enumerate(self.diffusion_models): if dm.name == name: - del self._diffusion_models[i] + del self.diffusion_models[i] self._generate_component_collections() return raise ValueError( f'No DiffusionModel with name {name} found. \n' - f'The available names are: {[dm.name for dm in self._diffusion_models]}' + f'The available names are: {[dm.name for dm in self.diffusion_models]}' ) def clear_diffusion_models(self) -> None: """Clear all DiffusionModels from the SampleModel.""" - self._diffusion_models.clear() + self.diffusion_models = [] self._generate_component_collections() # ------------------------------------------------------------------ @@ -214,6 +220,7 @@ def diffusion_models( self._diffusion_models = [] return if isinstance(value, DiffusionModelBase): + value.Q = self.Q self._diffusion_models = [value] return if not isinstance(value, list) or not all( @@ -223,6 +230,8 @@ def diffusion_models( 'diffusion_models must be a DiffusionModelBase, a list of DiffusionModelBase, ' 'or None' ) + for dm in value: + dm.Q = self.Q self._diffusion_models = value self._on_diffusion_models_change() @@ -492,8 +501,8 @@ def get_all_variables(self, Q_index: int | None = None) -> list[Parameter]: if self.temperature is not None: all_vars.append(self.temperature) - for diffusion_model in self._diffusion_models: - all_vars.extend(diffusion_model.get_all_variables()) + for diffusion_model in self.diffusion_models: + all_vars.extend(diffusion_model.get_all_variables(Q_index=Q_index)) return all_vars @@ -512,11 +521,8 @@ def _generate_component_collections(self) -> None: return # Generate components from diffusion models # and add to component collections - for diffusion_model in self._diffusion_models: - diffusion_collections = diffusion_model.create_component_collections( - Q=self.Q, - component_name=diffusion_model.name, - ) + for diffusion_model in self.diffusion_models: + diffusion_collections = diffusion_model.get_component_collections() for target, source in zip( self._component_collections, diffusion_collections, @@ -529,6 +535,15 @@ def _on_diffusion_models_change(self) -> None: """Handle changes to the diffusion models.""" self._generate_component_collections() + def _on_Q_change(self) -> None: + """Handle changes to the Q values.""" + + for diffusion_model in self.diffusion_models: + # This may be too aggressive + diffusion_model.clear_Q(confirm=True) + diffusion_model.Q = self.Q + self._generate_component_collections() + # ------------------------------------------------------------------ # dunder methods # ------------------------------------------------------------------ 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 4ffd2a5e..7941d406 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 @@ -43,7 +43,7 @@ def test_init_default(self, brownian_diffusion_model): ), ( { - 'unit': 123, + 'unit': 'meV', 'scale': 'invalid', 'diffusion_coefficient': 1.0, }, @@ -52,13 +52,31 @@ def test_init_default(self, brownian_diffusion_model): ), ( { - 'unit': 123, + 'unit': 'meV', + 'scale': -123.4, + 'diffusion_coefficient': 1.0, + }, + ValueError, + 'scale must be non-negative', + ), + ( + { + 'unit': 'meV', 'scale': 1.0, 'diffusion_coefficient': 'invalid', }, TypeError, 'diffusion_coefficient must be a number', ), + ( + { + 'unit': 'meV', + 'scale': 1.0, + 'diffusion_coefficient': -123.4, + }, + ValueError, + 'diffusion_coefficient must be non-negative', + ), ], ) def test_input_type_validation_raises(self, kwargs, expected_exception, expected_message): @@ -152,9 +170,9 @@ def test_calculate_QISF_type_error(self, brownian_diffusion_model): ) def test_create_component_collections(self, brownian_diffusion_model, Q): # WHEN - + brownian_diffusion_model.Q = Q # THEN - component_collections = brownian_diffusion_model.create_component_collections(Q=Q) + component_collections = brownian_diffusion_model.create_component_collections() # EXPECT expected_widths = brownian_diffusion_model.calculate_width(Q) @@ -166,36 +184,6 @@ def test_create_component_collections(self, brownian_diffusion_model, Q): assert np.isclose(component.width.value, expected_widths[model_index]) assert component.width.independent is False - def test_create_component_collections_component_name_must_be_string( - self, brownian_diffusion_model - ): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'component_name must be a string.'): - brownian_diffusion_model.create_component_collections( - Q=np.array([0.1, 0.2, 0.3]), component_name=123 - ) - - def test_create_component_collections_component_display_name_must_be_string( - self, brownian_diffusion_model - ): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'component_display_name must be a string.'): - brownian_diffusion_model.create_component_collections( - Q=np.array([0.1, 0.2, 0.3]), component_display_name=123 - ) - - def test_create_component_collections_Q_type_error(self, brownian_diffusion_model): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match='Q must be a '): - brownian_diffusion_model.create_component_collections(Q='invalid') # Invalid type - - def test_create_component_collections_Q_1dimensional_error(self, brownian_diffusion_model): - # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'Q must be a 1-dimensional array.'): - brownian_diffusion_model.create_component_collections( - Q=np.array([[0.1, 0.2], [0.3, 0.4]]) - ) # Invalid shape - 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 new file mode 100644 index 00000000..927b3442 --- /dev/null +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_delta_lorentz.py @@ -0,0 +1,422 @@ +# SPDX-FileCopyrightText: 2026 EasyScience contributors +# SPDX-License-Identifier: BSD-3-Clause + +import numpy as np +import pytest +from easyscience.variable import Parameter + +from easydynamics.sample_model.diffusion_model.delta_lorentz import DeltaLorentz + + +class TestDeltaLorentz: + @pytest.fixture + def delta_lorentz_model(self): + return DeltaLorentz() + + @pytest.fixture + def delta_lorentz_model_with_Q(self): + Q = np.linspace(0.5, 2, 7) + return DeltaLorentz( + Q=Q, + A_0=0.5, + lorentzian_width=0.0015, + allow_Q_variation={'A_0': True, 'lorentzian_width': True}, + ) + + def test_init_default(self, delta_lorentz_model): + # WHEN THEN EXPECT + assert delta_lorentz_model.display_name == 'DeltaLorentz' + assert delta_lorentz_model.unit == 'meV' + assert delta_lorentz_model.scale.value == pytest.approx(1.0) + assert delta_lorentz_model.mean_u_squared.value == pytest.approx(0.0) + assert delta_lorentz_model.A_0.value == pytest.approx(1.0) + assert delta_lorentz_model.lorentzian_width.value == pytest.approx(1.0) + + def test_init_with_Q(self, delta_lorentz_model_with_Q): + # WHEN THEN EXPECT + assert delta_lorentz_model_with_Q.display_name == 'DeltaLorentz' + assert delta_lorentz_model_with_Q.unit == 'meV' + assert delta_lorentz_model_with_Q.scale.value == pytest.approx(1.0) + assert delta_lorentz_model_with_Q.mean_u_squared.value == pytest.approx(0.0) + assert delta_lorentz_model_with_Q.A_0.value == pytest.approx(0.5) + assert delta_lorentz_model_with_Q.lorentzian_width.value == pytest.approx(0.0015) + assert delta_lorentz_model_with_Q._allow_Q_variation == { + 'A_0': True, + 'lorentzian_width': True, + } + assert len(delta_lorentz_model_with_Q._A_0_list) == len(delta_lorentz_model_with_Q.Q) + assert len(delta_lorentz_model_with_Q._lorentzian_width_list) == len( + delta_lorentz_model_with_Q.Q + ) + assert all(pytest.approx(a.value) == 0.5 for a in delta_lorentz_model_with_Q._A_0_list) + assert all( + pytest.approx(lw.value) == 0.0015 + for lw in delta_lorentz_model_with_Q._lorentzian_width_list + ) + + @pytest.mark.parametrize( + 'kwargs,expected_exception, expected_message', + [ + ( + { + 'mean_u_squared': -1.0, + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + ValueError, + 'mean_u_squared must be non-negative', + ), + ( + { + 'mean_u_squared': 'not a number', + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'mean_u_squared must be a number', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': -1.0, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + ValueError, + 'A_0 must be between 0 and 1', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 'not a number', + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'A_0 must be a number', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': -1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + ValueError, + 'lorentzian_width must be ', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': 'not a number', + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'lorentzian_width must be a number', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': 'Not a dict', + 'delta_name': 'Delta', + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'allow_Q_variation must be a dict', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 123, + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'delta_name must be a string', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': None, + 'delta_display_name': 'DeltaDisplay', + }, + TypeError, + 'delta_name must be a string', + ), + ( + { + 'mean_u_squared': 0.1, + 'A_0': 0.5, + 'lorentzian_width': 1.0, + 'allow_Q_variation': {'A_0': True, 'lorentzian_width': True}, + 'delta_name': 'Delta', + 'delta_display_name': 123, + }, + TypeError, + 'delta_display_name must be a string', + ), + ], + ids=[ + 'mean_u_squared negative', + 'mean_u_squared not a number', + 'A_0 negative', + 'A_0 not a number', + 'lorentzian_width negative', + 'lorentzian_width not a number', + 'allow_Q_variation not a dict', + 'delta_name not a string', + 'delta_name not a string (None)', + 'delta_display_name not a string', + ], + ) + def test_input_type_validation_raises(self, kwargs, expected_exception, expected_message): + with pytest.raises(expected_exception, match=expected_message): + DeltaLorentz(**kwargs) + + # ------------------------------------------------------------------ + # Properties + # ------------------------------------------------------------------ + @pytest.mark.parametrize( + ('attribute', 'value', 'expected'), + [ + ('mean_u_squared', 2.0, 2.0), + ('mean_u_squared', 0.0, 0.0), + ('mean_u_squared', 5, 5.0), + ('A_0', 0.0, 0.0), + ('A_0', 1.0, 1.0), + ('A_0', 0.5, 0.5), + ('lorentzian_width', 1.5, 1.5), + ('delta_name', 'delta', 'delta'), + ('delta_display_name', 'display', 'display'), + ('delta_display_name', None, None), + ], + ids=[ + 'mean_u_squared set to 2.0', + 'mean_u_squared set to 0.0', + 'mean_u_squared set to 5 (int)', + 'A_0 set to 0.0', + 'A_0 set to 1.0', + 'A_0 set to 0.5', + 'lorentzian_width set to 1.5', + "delta_name set to 'delta'", + "delta_display_name set to 'display'", + 'delta_display_name set to None', + ], + ) + def test_setters_valid( + self, + delta_lorentz_model, + attribute, + value, + expected, + ): + # WHEN + setattr(delta_lorentz_model, attribute, value) + + # THEN + result = getattr(delta_lorentz_model, attribute) + + # Handle Parameters + if isinstance(result, Parameter): + result = result.value + + # EXPECT + assert result == expected + + @pytest.mark.parametrize( + ('attribute', 'value', 'exception', 'message'), + [ + ( + 'mean_u_squared', + -1.0, + ValueError, + r'mean_u_squared must be non-negative.', + ), + ( + 'mean_u_squared', + 'invalid', + TypeError, + r'mean_u_squared must be a number.', + ), + ( + 'A_0', + -0.1, + ValueError, + r'A_0 must be between 0 and 1.', + ), + ( + 'A_0', + 1.1, + ValueError, + r'A_0 must be between 0 and 1.', + ), + ( + 'A_0', + 'invalid', + TypeError, + r'A_0 must be a number.', + ), + ( + 'A_1', + 0.5, + AttributeError, + r'A_1 is a dependent parameter and cannot be set directly.', + ), + ( + 'lorentzian_width', + -0.1, + ValueError, + r'lorentzian_width must be.', + ), + ( + 'lorentzian_width', + 'invalid', + TypeError, + r'lorentzian_width must be a number.', + ), + ( + 'delta_name', + 1, + TypeError, + r'delta_name must be a string.', + ), + ( + 'delta_name', + None, + TypeError, + r'delta_name must be a string.', + ), + ( + 'delta_display_name', + 1, + TypeError, + r'delta_display_name must be a string or None.', + ), + ( + 'delta_display_name', + [], + TypeError, + r'delta_display_name must be a string or None.', + ), + ], + ids=[ + 'mean_u_squared negative', + 'mean_u_squared not a number', + 'A_0 less than 0', + 'A_0 greater than 1', + 'A_0 not a number', + 'A_1 set directly', + 'lorentzian_width negative', + 'lorentzian_width not a number', + 'delta_name not a string', + 'delta_name not a string (None)', + 'delta_display_name not a string', + 'delta_display_name not a string (list)', + ], + ) + def test_setters_invalid( + self, + delta_lorentz_model, + attribute, + value, + exception, + message, + ): + # WHEN THEN EXPECT + 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 + width = delta_lorentz_model.calculate_width(Q=0.5) + + # EXPECT + assert len(width) == 1 + assert width[0] == pytest.approx(1.0) + + def test_calculate_width_with_Q(self, delta_lorentz_model_with_Q): + # WHEN THEN + width = delta_lorentz_model_with_Q.calculate_width() + + # EXPECT + assert len(width) == len(delta_lorentz_model_with_Q.Q) + assert all(width_i == pytest.approx(0.0015) for width_i in width) + + def test_calculate_EISF(self, delta_lorentz_model): + # WHEN + + # THEN + eisf = delta_lorentz_model.calculate_EISF(Q=0.5) + + # EXPECT + assert len(eisf) == 1 + expected = delta_lorentz_model.A_0.value * np.exp( + -delta_lorentz_model.mean_u_squared.value * 0.5**2 + ) + assert eisf[0] == pytest.approx(expected) + + def test_calculate_EISF_with_Q(self, delta_lorentz_model_with_Q): + # WHEN + + # THEN + eisf = delta_lorentz_model_with_Q.calculate_EISF() + + # EXPECT + assert len(eisf) == len(delta_lorentz_model_with_Q.Q) + for i in range(len(eisf)): + expected = delta_lorentz_model_with_Q._A_0_list[i].value * np.exp( + -delta_lorentz_model_with_Q.mean_u_squared.value + * delta_lorentz_model_with_Q.Q[i] ** 2 + ) + assert eisf[i] == pytest.approx(expected) + + def test_calculate_QISF(self, delta_lorentz_model): + # WHEN THEN + qisf = delta_lorentz_model.calculate_QISF(Q=0.5) + + # EXPECT + assert len(qisf) == 1 + expected = delta_lorentz_model.A_1.value * np.exp( + -delta_lorentz_model.mean_u_squared.value * 0.5**2 + ) + assert qisf[0] == pytest.approx(expected) + + def test_calculate_QISF_with_Q(self, delta_lorentz_model_with_Q): + # WHEN THEN + qisf = delta_lorentz_model_with_Q.calculate_QISF() + + # EXPECT + assert len(qisf) == len(delta_lorentz_model_with_Q.Q) + for i in range(len(qisf)): + expected = delta_lorentz_model_with_Q._A_1_list[i].value * np.exp( + -delta_lorentz_model_with_Q.mean_u_squared.value + * delta_lorentz_model_with_Q.Q[i] ** 2 + ) + + assert qisf[i] == pytest.approx(expected) diff --git a/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model.py b/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model.py index c46dc561..2caf0b90 100644 --- a/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model.py +++ b/tests/unit/easydynamics/sample_model/diffusion_model/test_diffusion_model.py @@ -1,7 +1,9 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +import numpy as np import pytest +from easyscience.variable.parameter import Parameter from easydynamics.sample_model.diffusion_model.diffusion_model_base import DiffusionModelBase @@ -9,13 +11,24 @@ class TestDiffusionModel: @pytest.fixture def diffusion_model(self): - return DiffusionModelBase(display_name='TestDiffusionModel', unit='meV') + return DiffusionModelBase() def test_init_default(self, diffusion_model): # WHEN THEN EXPECT - assert diffusion_model.display_name == 'TestDiffusionModel' + assert diffusion_model.display_name == 'DiffusionModel' + assert diffusion_model.name == 'DiffusionModel' + assert diffusion_model.lorentzian_name == 'DiffusionModel' + assert diffusion_model.lorentzian_display_name == 'DiffusionModel' assert diffusion_model.unit == 'meV' + def test_init_raises(self): + # WHEN THEN EXPECT + with pytest.raises(TypeError, match=r'lorentzian_name must be a string'): + DiffusionModelBase(lorentzian_name=123) + + with pytest.raises(TypeError, match=r'lorentzian_display_name must be a string or None'): + DiffusionModelBase(lorentzian_display_name=123) + def test_unit_setter_raises(self, diffusion_model): # WHEN THEN EXPECT with pytest.raises( @@ -24,22 +37,115 @@ def test_unit_setter_raises(self, diffusion_model): ): diffusion_model.unit = 'eV' - def test_scale_setter(self, diffusion_model): + @pytest.mark.parametrize( + ('attribute', 'value', 'expected'), + [ + ('scale', 2.0, 2.0), + ('scale', 0.0, 0.0), + ('scale', 5, 5.0), + ('lorentzian_name', 'lorentzian', 'lorentzian'), + ('lorentzian_name', '', ''), + ('lorentzian_display_name', 'display', 'display'), + ('lorentzian_display_name', None, None), + ], + ) + def test_setters_valid( + self, + diffusion_model, + attribute, + value, + expected, + ): # WHEN - diffusion_model.scale = 2.0 - # THEN EXPECT - assert diffusion_model.scale.value == pytest.approx(2.0) + # THEN + setattr(diffusion_model, attribute, value) - def test_scale_setter_negative_raises(self, diffusion_model): + # EXPECT + result = getattr(diffusion_model, attribute) + + # Handle Parameters + if isinstance(result, Parameter): + result = result.value + + assert result == expected + + @pytest.mark.parametrize( + ('attribute', 'value', 'exception', 'message'), + [ + ( + 'scale', + -1.0, + ValueError, + r'scale must be non-negative.', + ), + ( + 'scale', + 'invalid', + TypeError, + r'scale must be a number.', + ), + ( + 'lorentzian_name', + 1, + TypeError, + r'lorentzian_name must be a string.', + ), + ( + 'lorentzian_name', + None, + TypeError, + r'lorentzian_name must be a string.', + ), + ( + 'lorentzian_display_name', + 1, + TypeError, + r'lorentzian_display_name must be a string or None.', + ), + ( + 'lorentzian_display_name', + [], + TypeError, + r'lorentzian_display_name must be a string or None.', + ), + ], + ) + def test_setters_invalid( + self, + diffusion_model, + attribute, + value, + exception, + message, + ): # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'scale must be non-negative.'): - diffusion_model.scale = -1.0 # Invalid negative value + with pytest.raises(exception, match=message): + setattr(diffusion_model, attribute, value) - def test_scale_setter_raises(self, diffusion_model): + def test_Q_property(self, diffusion_model): # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'scale must be a number.'): - diffusion_model.scale = 'invalid' # Invalid type + assert diffusion_model.Q is None + + # THEN + diffusion_model.Q = [1.0, 2.0, 3.0] + + # EXPECT + np.testing.assert_allclose(diffusion_model.Q, [1.0, 2.0, 3.0]) + + # THEN EXPECT + with pytest.raises(ValueError, match=r'New Q values are not similar to the old ones'): + diffusion_model.Q = [10.0, 20.0, 30.0] + + # THEN EXPECT + with pytest.raises(ValueError, match=r'Clearing Q values requires confirmation'): + diffusion_model.clear_Q() + + # THEN + diffusion_model.clear_Q(confirm=True) + + # EXPECT + assert diffusion_model.Q is None def test_repr(self, diffusion_model): # WHEN THEN @@ -47,5 +153,4 @@ def test_repr(self, diffusion_model): # EXPECT assert 'DiffusionModelBase' in repr_str - assert 'display_name=TestDiffusionModel' in repr_str assert 'unit=meV' in repr_str 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 d1fbf8ff..014a221f 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 @@ -63,6 +63,16 @@ def test_init_default(self, jump_diffusion_model): TypeError, 'diffusion_coefficient must be a number', ), + ( + { + 'unit': 'meV', + 'scale': 1.0, + 'diffusion_coefficient': -1.0, + 'relaxation_time': 1.0, + }, + ValueError, + 'diffusion_coefficient must be non-negative', + ), ( { 'unit': 'meV', @@ -73,6 +83,16 @@ def test_init_default(self, jump_diffusion_model): TypeError, 'relaxation_time must be a number', ), + ( + { + 'unit': 'meV', + 'scale': 1.0, + 'diffusion_coefficient': 1.0, + 'relaxation_time': -1.0, + }, + ValueError, + 'relaxation_time must be non-negative', + ), ], ) def test_input_type_validation_raises(self, kwargs, expected_exception, expected_message): @@ -190,9 +210,10 @@ def test_calculate_QISF_type_error(self, jump_diffusion_model): ) def test_create_component_collections(self, jump_diffusion_model, Q): # WHEN + jump_diffusion_model.Q = Q # THEN - component_collections = jump_diffusion_model.create_component_collections(Q=Q) + component_collections = jump_diffusion_model.create_component_collections() # EXPECT expected_widths = jump_diffusion_model.calculate_width(Q) @@ -204,36 +225,6 @@ def test_create_component_collections(self, jump_diffusion_model, Q): assert np.isclose(component.width.value, expected_widths[model_index]) assert component.width.independent is False - def test_create_component_collections_component_name_must_be_string( - self, jump_diffusion_model - ): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'component_name must be a string.'): - jump_diffusion_model.create_component_collections( - Q=np.array([0.1, 0.2, 0.3]), component_name=123 - ) - - def test_create_component_collections_component_display_name_must_be_string( - self, jump_diffusion_model - ): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match=r'component_display_name must be a string.'): - jump_diffusion_model.create_component_collections( - Q=np.array([0.1, 0.2, 0.3]), component_display_name=123 - ) - - def test_create_component_collections_Q_type_error(self, jump_diffusion_model): - # WHEN THEN EXPECT - with pytest.raises(TypeError, match='Q must be a '): - jump_diffusion_model.create_component_collections(Q='invalid') # Invalid type - - def test_create_component_collections_Q_1dimensional_error(self, jump_diffusion_model): - # WHEN THEN EXPECT - with pytest.raises(ValueError, match=r'Q must be a 1-dimensional array.'): - jump_diffusion_model.create_component_collections( - Q=np.array([[0.1, 0.2], [0.3, 0.4]]) - ) # Invalid shape - 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_sample_model.py b/tests/unit/easydynamics/sample_model/test_sample_model.py index 337266f5..3055440c 100644 --- a/tests/unit/easydynamics/sample_model/test_sample_model.py +++ b/tests/unit/easydynamics/sample_model/test_sample_model.py @@ -469,7 +469,6 @@ def test_generate_component_collections(self, sample_model): assert collection[0].area.value == pytest.approx(1.0) assert collection[1].name == 'TestLorentzian1Name' assert collection[1].area.value == pytest.approx(2.0) - assert collection[2].name == 'DiffusionModelName' assert isinstance(collection[2], Lorentzian) def test_get_all_variables(self, sample_model):