diff --git a/docs/docs/tutorials/sample_model.ipynb b/docs/docs/tutorials/sample_model.ipynb index 882048ce..002bbfca 100644 --- a/docs/docs/tutorials/sample_model.ipynb +++ b/docs/docs/tutorials/sample_model.ipynb @@ -148,7 +148,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.4" + "version": "3.14.5" } }, "nbformat": 4, diff --git a/docs/docs/tutorials/tutorial0_more_advanced.ipynb b/docs/docs/tutorials/tutorial0_more_advanced.ipynb index 25088780..11e6e7d5 100644 --- a/docs/docs/tutorials/tutorial0_more_advanced.ipynb +++ b/docs/docs/tutorials/tutorial0_more_advanced.ipynb @@ -366,7 +366,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.4" + "version": "3.14.5" } }, "nbformat": 4, diff --git a/docs/docs/tutorials/tutorial1_brownian.ipynb b/docs/docs/tutorials/tutorial1_brownian.ipynb index 40a6c21f..269bb259 100644 --- a/docs/docs/tutorials/tutorial1_brownian.ipynb +++ b/docs/docs/tutorials/tutorial1_brownian.ipynb @@ -94,20 +94,7 @@ "id": "6c87b01c", "metadata": {}, "source": [ - "We now want to fit the vanadium data to determine our resolution. The scattering from vanadium is almost exclusively incoherent elastic, so we model it as a delta function. We do this by creating a `SampleModel` and adding a `DeltaFunction` component to it. The component acts as a template and gets copied to every `Q` when we attach the `SampleModel` to our `Analysis` object. Let's create the `SampleModel`.\n", - "\n", - "We do not give the `DeltaFunction` a `center` value. In this case, the center will be fixed at 0 energy transfer. We set the start value of the area to 1." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6762faba", - "metadata": {}, - "outputs": [], - "source": [ - "delta_function = sm.DeltaFunction(name='DeltaFunction', area=1)\n", - "sample_model = sm.SampleModel(components=delta_function)" + "We now want to fit the vanadium data to determine our resolution. The scattering from vanadium is almost exclusively incoherent elastic, so it can be modelled as a delta function. There are two ways to do this. The first is to add a `DeltaFunction` component to a `SampleModel` and next create a `ResolutionModel`. The second method, which we here shall use, is to fit the data to a `SampleModel` and then convert it to a `ResolutionModel`." ] }, { @@ -115,11 +102,7 @@ "id": "dc82774e", "metadata": {}, "source": [ - "We now want to define our resolution function. We will here model it as a Gaussian. We create a `ComponentCollection` and append the `Gaussian` to it. We can add as many components to our resolution as we like; sometimes you need several Gaussians and other functions to accurately describe the resolution.\n", - "\n", - "We fix the area of the resolution to have value 1. If we did not do this, we would fit both the area of the delta function and of the resolution Gaussian, and the fit would never converge.\n", - "\n", - "We finally insert the components in a `ResolutionModel`" + "We now want to define our resolution function. We will here model it as a Gaussian. We create a `ComponentCollection` and append the `Gaussian` to it. We can add as many components to our resolution as we like; sometimes you need several Gaussians and other functions to accurately describe the resolution." ] }, { @@ -129,11 +112,10 @@ "metadata": {}, "outputs": [], "source": [ - "resolution_components = sm.ComponentCollection()\n", + "vanadium_components = sm.ComponentCollection()\n", "res_gauss = sm.Gaussian(width=0.1, area=1, name='Res. Gauss')\n", - "res_gauss.area.fixed = True\n", - "resolution_components.append_component(res_gauss)\n", - "resolution_model = sm.ResolutionModel(components=resolution_components)" + "vanadium_components.append_component(res_gauss)\n", + "vanadium_model = sm.SampleModel(components=vanadium_components)" ] }, { @@ -141,7 +123,7 @@ "id": "088ac17d", "metadata": {}, "source": [ - "The background intensity was not 0, so we also create a background model. We use a `Polynomial` with a single coefficient, i.e. a flat background. We here show how to create the `BackgroundModel` and add the background in a single line. We could of course also add it like we did for the `SampleModel` or first create a `ComponentCollection` like we did for the `ResolutionModel`" + "The background intensity was not 0, so we also create a background model. We use a `Polynomial` with a single coefficient, i.e. a flat background. We here show how to create the `BackgroundModel` and add the background in a single line. We could of course also add it like we did for the `SampleModel` or first create a `ComponentCollection` like we just did for the resolution. " ] }, { @@ -159,7 +141,7 @@ "id": "eae3d14b", "metadata": {}, "source": [ - "We combine the resolution abd background model into an `InstrumentModel`. This model also contains a fittable energy offset to account for instrument misalignment. All components are centered at this energy offset." + "We add the background model to an `InstrumentModel`. This model also contains a fittable energy offset to account for instrument misalignment. All components are centered at this energy offset." ] }, { @@ -170,7 +152,6 @@ "outputs": [], "source": [ "instrument_model = sm.InstrumentModel(\n", - " resolution_model=resolution_model,\n", " background_model=background_model,\n", ")" ] @@ -193,7 +174,7 @@ "vanadium_analysis = edyn.Analysis(\n", " display_name='Vanadium Full Analysis',\n", " experiment=vanadium_experiment,\n", - " sample_model=sample_model,\n", + " sample_model=vanadium_model,\n", " instrument_model=instrument_model,\n", ")" ] @@ -265,7 +246,7 @@ "outputs": [], "source": [ "# Plot some of fitted parameters as a function of Q\n", - "vanadium_analysis.plot_parameters(names=['DeltaFunction area'])" + "vanadium_analysis.plot_parameters(names=['Res. Gauss area'])" ] }, { @@ -356,7 +337,7 @@ "id": "927b8fb5", "metadata": {}, "source": [ - "We also create a new instrument_model and attach it to our analysis, giving it the resolution model determined in the vanadium analysis. We further fix all parameters in the resolution model and normalize it." + "We also create a new instrument_model and attach it to our analysis. We now give it the sample model from the vanadium fit. All parameters are automatically fixed and the resolution model is normalized to have area 1. " ] }, { @@ -368,10 +349,8 @@ "source": [ "instrument_model = sm.InstrumentModel(\n", " background_model=background_model,\n", - " resolution_model=vanadium_analysis.instrument_model.resolution_model,\n", + " resolution_model=vanadium_analysis.sample_model,\n", ")\n", - "instrument_model.resolution_model.fix_all_parameters()\n", - "instrument_model.normalize_resolution()\n", "\n", "diffusion_analysis = edyn.Analysis(\n", " display_name='Diffusion Analysis',\n", @@ -607,7 +586,7 @@ "source": [ "instrument_model = sm.InstrumentModel(\n", " background_model=background_model,\n", - " resolution_model=vanadium_analysis.instrument_model.resolution_model,\n", + " resolution_model=vanadium_analysis.sample_model,\n", ")" ] }, @@ -726,7 +705,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.4" + "version": "3.14.5" } }, "nbformat": 4, diff --git a/docs/docs/tutorials/tutorial2_nanoparticles.ipynb b/docs/docs/tutorials/tutorial2_nanoparticles.ipynb index 2fb4fc4a..22de2490 100644 --- a/docs/docs/tutorials/tutorial2_nanoparticles.ipynb +++ b/docs/docs/tutorials/tutorial2_nanoparticles.ipynb @@ -123,7 +123,7 @@ "source": [ "The resolution function can to a good approximation be modelled as a Gaussian. In truth, it has small Lorentzian tails, but we leave it as an exercise for the interested reader to add this.\n", "\n", - "We define the `SampleModel` to be a `DeltaFunction`, since there is essentially only elastic scattering present. We also add a constant background using a `Polynomial`. \n", + "We fit the resolution directly as in the previous tutorial, since there is essentially only elastic scattering present. We also add a constant background using a `Polynomial`. \n", "\n", "As in Tutorial 1 we place everything in an `Analysis` object and plot the start guesses." ] @@ -135,16 +135,12 @@ "metadata": {}, "outputs": [], "source": [ - "delta_function = sm.DeltaFunction(area=100)\n", - "res_sample_model = sm.SampleModel(components=delta_function)\n", - "\n", - "res_resolution_model = sm.ResolutionModel()\n", + "res_sample_model = sm.SampleModel()\n", "res_components = sm.ComponentCollection()\n", - "res_gauss = sm.Gaussian(area=1, width=0.02)\n", - "res_gauss.area.fixed = True\n", + "res_gauss = sm.Gaussian(area=40, width=0.02)\n", "\n", "res_components.append_component(res_gauss)\n", - "res_resolution_model.components = res_components\n", + "res_sample_model.components = res_components\n", "\n", "background_model = sm.BackgroundModel()\n", "polynomial = sm.Polynomial(coefficients=[1.5])\n", @@ -153,7 +149,6 @@ "\n", "\n", "res_instrument_model = sm.InstrumentModel(\n", - " resolution_model=res_resolution_model,\n", " background_model=background_model,\n", ")\n", "\n", @@ -267,10 +262,8 @@ "\n", "instrument_model = sm.InstrumentModel(\n", " background_model=background_model,\n", - " resolution_model=res_analysis.instrument_model.resolution_model,\n", + " resolution_model=res_analysis.sample_model,\n", ")\n", - "instrument_model.resolution_model.fix_all_parameters()\n", - "instrument_model.normalize_resolution()\n", "\n", "\n", "analysis = edyn.Analysis(\n", @@ -407,10 +400,8 @@ "\n", "instrument_model = sm.InstrumentModel(\n", " background_model=background_model,\n", - " resolution_model=res_analysis.instrument_model.resolution_model,\n", + " resolution_model=res_analysis.sample_model,\n", ")\n", - "instrument_model.resolution_model.fix_all_parameters()\n", - "instrument_model.normalize_resolution()\n", "\n", "# Create the analysis object\n", "mag_analysis = edyn.Analysis(\n", @@ -616,7 +607,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.4" + "version": "3.14.5" } }, "nbformat": 4, diff --git a/src/easydynamics/sample_model/instrument_model.py b/src/easydynamics/sample_model/instrument_model.py index 6468cd3b..fd0ae458 100644 --- a/src/easydynamics/sample_model/instrument_model.py +++ b/src/easydynamics/sample_model/instrument_model.py @@ -10,6 +10,7 @@ from easydynamics.sample_model.background_model import BackgroundModel from easydynamics.sample_model.resolution_model import ResolutionModel +from easydynamics.sample_model.sample_model import SampleModel from easydynamics.utils.utils import Numeric from easydynamics.utils.utils import Q_type from easydynamics.utils.utils import _validate_and_convert_Q @@ -29,7 +30,7 @@ def __init__( display_name: str = 'MyInstrumentModel', unique_name: str | None = None, Q: Q_type | None = None, - resolution_model: ResolutionModel | None = None, + resolution_model: ResolutionModel | SampleModel | None = None, background_model: BackgroundModel | None = None, energy_offset: Numeric | None = None, unit: str | sc.Unit = 'meV', @@ -45,9 +46,10 @@ def __init__( The unique name of the InstrumentModel. Q : Q_type | None, default=None The Q values where the instrument is modelled. - resolution_model : ResolutionModel | None, default=None - The resolution model of the instrument. If None, an empty resolution model is created - and no resolution convolution is carried out. + resolution_model : ResolutionModel | SampleModel | None, default=None + The resolution model of the instrument. If a SampleModel it will be converted to a + ResolutionModel. If None, an empty resolution model is created and no resolution + convolution is carried out. background_model : BackgroundModel | None, default=None The background model of the instrument. If None, an empty background model is created, and the background evaluates to 0. @@ -73,11 +75,13 @@ def __init__( if resolution_model is None: self._resolution_model = ResolutionModel() else: - if not isinstance(resolution_model, ResolutionModel): + if not isinstance(resolution_model, (ResolutionModel, SampleModel)): raise TypeError( - f'resolution_model must be a ResolutionModel or None, ' + f'resolution_model must be a ResolutionModel, a SampleModel or None, ' f'got {type(resolution_model).__name__}' ) + if isinstance(resolution_model, SampleModel): + resolution_model = ResolutionModel.from_sample_model(resolution_model) self._resolution_model = resolution_model if background_model is None: @@ -122,24 +126,27 @@ def resolution_model(self) -> ResolutionModel: return self._resolution_model @resolution_model.setter - def resolution_model(self, value: ResolutionModel) -> None: + def resolution_model(self, value: ResolutionModel | SampleModel) -> None: """ Set the resolution model of the instrument. Parameters ---------- - value : ResolutionModel + value : ResolutionModel | SampleModel The new resolution model of the instrument. Raises ------ TypeError - If value is not a ResolutionModel. + If value is not a ResolutionModel or SampleModel. """ - if not isinstance(value, ResolutionModel): + if not isinstance(value, (ResolutionModel, SampleModel)): raise TypeError( - f'resolution_model must be a ResolutionModel, got {type(value).__name__}' + f'resolution_model must be a ResolutionModel or SampleModel, ' + f'got {type(value).__name__}' ) + if isinstance(value, SampleModel): + value = ResolutionModel.from_sample_model(value) self._resolution_model = value self._on_resolution_model_change() diff --git a/src/easydynamics/sample_model/resolution_model.py b/src/easydynamics/sample_model/resolution_model.py index 4f35cc30..b90f2b54 100644 --- a/src/easydynamics/sample_model/resolution_model.py +++ b/src/easydynamics/sample_model/resolution_model.py @@ -1,13 +1,17 @@ # SPDX-FileCopyrightText: 2026 EasyScience contributors # SPDX-License-Identifier: BSD-3-Clause +from copy import copy + import scipp as sc from easydynamics.sample_model.component_collection import ComponentCollection from easydynamics.sample_model.components import DeltaFunction from easydynamics.sample_model.components import Polynomial +from easydynamics.sample_model.components.exponential import Exponential from easydynamics.sample_model.components.model_component import ModelComponent from easydynamics.sample_model.model_base import ModelBase +from easydynamics.sample_model.sample_model import SampleModel from easydynamics.utils.utils import Q_type @@ -29,11 +33,11 @@ def __init__( Parameters ---------- - display_name : str, default='MyResolutionModel' + display_name : str, default="MyResolutionModel" Display name of the model. unique_name : str | None, default=None Unique name of the model. If None, a unique name will be generated. - unit : str | sc.Unit, default='meV' + unit : str | sc.Unit, default="meV" Unit of the model. components : ModelComponent | ComponentCollection | None, default=None Template components of the model. If None, no components are added. These components @@ -70,9 +74,70 @@ def append_component(self, component: ModelComponent | ComponentCollection) -> N components = component if isinstance(component, ComponentCollection) else (component,) for comp in components: - if isinstance(comp, (DeltaFunction, Polynomial)): + if isinstance(comp, (DeltaFunction, Polynomial, Exponential)): raise TypeError( f'Component in ResolutionModel cannot be a {comp.__class__.__name__}' ) super().append_component(component) + + @classmethod + def from_sample_model( + cls, + sample_model: SampleModel, + normalize_area: bool = True, + fix_parameters: bool = True, + ) -> 'ResolutionModel': + """ + Create a ResolutionModel from a SampleModel. + + Parameters + ---------- + sample_model : SampleModel + SampleModel to create the ResolutionModel from. + normalize_area : bool, default=True + Whether to normalize the components in the ResolutionModel to have area 1. + fix_parameters : bool, default=True + Whether to fix the parameters in the ResolutionModel. + + Returns + ------- + "ResolutionModel" + ResolutionModel created from the SampleModel. + + Raises + ------ + TypeError + If sample_model is not a SampleModel, or if normalize_area or fix_parameters are not + bool. + """ + if not isinstance(sample_model, SampleModel): + raise TypeError( + f'sample_model must be an instance of SampleModel. Got {type(sample_model)}.' + ) + + if not isinstance(normalize_area, bool): + raise TypeError('normalize_area must be True or False.') + + if not isinstance(fix_parameters, bool): + raise TypeError('fix_parameters must be True or False.') + + resolution_model = cls( + display_name=sample_model.display_name, + unit=sample_model.unit, + components=sample_model.components, + Q=sample_model.Q, + ) + + if sample_model.Q is not None: + for index in range(len(sample_model.Q)): + resolution_model._component_collections[index] = copy( + sample_model.get_component_collection(Q_index=index) + ) + if normalize_area: + resolution_model.normalize_area() + + if fix_parameters: + resolution_model.fix_all_parameters() + + return resolution_model diff --git a/tests/unit/easydynamics/sample_model/test_instrument_model.py b/tests/unit/easydynamics/sample_model/test_instrument_model.py index 781d0260..a61acff1 100644 --- a/tests/unit/easydynamics/sample_model/test_instrument_model.py +++ b/tests/unit/easydynamics/sample_model/test_instrument_model.py @@ -13,6 +13,7 @@ from easydynamics.sample_model.background_model import BackgroundModel from easydynamics.sample_model.instrument_model import InstrumentModel from easydynamics.sample_model.resolution_model import ResolutionModel +from easydynamics.sample_model.sample_model import SampleModel class TestInstrumentModel: @@ -56,6 +57,11 @@ def background_model(self): component = Polynomial(coefficients=[1.0, 2.0]) return BackgroundModel(components=component) + @pytest.fixture + def sample_model(self): + component = Gaussian() + return SampleModel(components=component) + def test_init(self, instrument_model): # WHEN THEN model = instrument_model @@ -79,6 +85,15 @@ def test_init_defaults(self): assert isinstance(model.resolution_model, ResolutionModel) assert model.Q is None + def test_init_sample_model_as_resolution_model(self, sample_model): + # WHEN THEN + model = InstrumentModel(resolution_model=sample_model) + + # EXPECT + assert model.display_name == 'MyInstrumentModel' + assert isinstance(model.background_model, BackgroundModel) + assert isinstance(model.resolution_model, ResolutionModel) + @pytest.mark.parametrize( 'kwargs, expected_exception, expected_message', [ @@ -127,6 +142,32 @@ def test_resolution_model_setter_calls_update(self, instrument_model, resolution assert instrument_model._resolution_model is resolution_model instrument_model._on_resolution_model_change.assert_called_once() + def test_resolution_model_setter_converts_sample_model( + self, + instrument_model, + ): + """ + Test that setting a SampleModel as the resolution_model calls the method to convert it to + a ResolutionModel, and updates the resolution_model and calls the change callback. + """ + # WHEN + sample_model = SampleModel() + converted_model = ResolutionModel() + instrument_model._on_resolution_model_change = MagicMock() + + with patch.object( + ResolutionModel, + 'from_sample_model', + return_value=converted_model, + ) as mock_from_sample_model: + # THEN + instrument_model.resolution_model = sample_model + + # EXPECT + mock_from_sample_model.assert_called_once_with(sample_model) + assert instrument_model._resolution_model is converted_model + instrument_model._on_resolution_model_change.assert_called_once_with() + def test_resolution_model_setter_raises(self, instrument_model): # WHEN / THEN / EXPECT with pytest.raises( diff --git a/tests/unit/easydynamics/sample_model/test_resolution_model.py b/tests/unit/easydynamics/sample_model/test_resolution_model.py index 9febca15..f894daf2 100644 --- a/tests/unit/easydynamics/sample_model/test_resolution_model.py +++ b/tests/unit/easydynamics/sample_model/test_resolution_model.py @@ -10,6 +10,7 @@ from easydynamics.sample_model import Lorentzian from easydynamics.sample_model import Polynomial from easydynamics.sample_model.resolution_model import ResolutionModel +from easydynamics.sample_model.sample_model import SampleModel class TestResolutionModel: @@ -39,6 +40,36 @@ def resolution_model(self): Q=np.array([1.0, 2.0, 3.0]), ) + @pytest.fixture + def sample_model(self): + component1 = Gaussian( + name='TestGaussian1Name', + display_name='TestGaussian1Display', + area=1.0, + center=0.0, + width=1.0, + unit='meV', + ) + component2 = Lorentzian( + name='TestLorentzian1Name', + display_name='TestLorentzian1Display', + area=2.0, + center=1.0, + width=0.5, + unit='meV', + ) + component_collection = ComponentCollection() + component_collection.append_component(component1) + component_collection.append_component(component2) + + return SampleModel( + display_name='InitModel', + components=component_collection, + unit='meV', + Q=np.array([1.0, 2.0, 3.0]), + temperature=10.0, + ) + def test_init(self, resolution_model): # WHEN THEN model = resolution_model @@ -151,3 +182,122 @@ def test_append_invalid_component_type_raises(self, resolution_model, invalid_co collection = ComponentCollection() collection.append_component(invalid_component) resolution_model.append_component(collection) + + @pytest.mark.parametrize( + ('kwargs', 'expected_area_sum', 'all_fixed'), + [ + ({}, 1.0, True), + ( + {'normalize_area': False, 'fix_parameters': False}, + 3.0, + False, + ), + ], + ids=[ + 'default_kwargs', + 'no_normalization_no_fixing', + ], + ) + def test_from_sample_model( + self, + sample_model, + kwargs, + expected_area_sum, + all_fixed, + ): + # WHEN + + # THEN + resolution_model = ResolutionModel.from_sample_model( + sample_model=sample_model, + **kwargs, + ) + + # EXPECT + assert resolution_model.display_name == 'InitModel' + assert resolution_model.unit == 'meV' + assert len(resolution_model.components) == 2 + np.testing.assert_array_equal( + resolution_model.Q, + np.array([1.0, 2.0, 3.0]), + ) + + for Q_index in range(len(resolution_model.Q)): + assert resolution_model.get_component_collection( + Q_index + ) is not sample_model.get_component_collection(Q_index) + + res_components = resolution_model.get_component_collection(Q_index) + assert (res_components[0].area.value + res_components[1].area.value) == pytest.approx( + expected_area_sum + ) + + sample_components = sample_model.get_component_collection(Q_index) + for res_comp, sample_comp in zip(res_components, sample_components, strict=True): + assert res_comp.name == sample_comp.name + assert res_comp.display_name == sample_comp.display_name + assert res_comp.center.value == sample_comp.center.value + assert res_comp.width.value == sample_comp.width.value + + variables = resolution_model.get_all_variables() + assert all(var.fixed for var in variables) is all_fixed + + def test_from_sample_model_with_no_Q(self, sample_model): + # WHEN + sample_model_no_Q = SampleModel( + components=sample_model.components, + ) + resolution_model = ResolutionModel.from_sample_model(sample_model_no_Q) + + # THEN / EXPECT + assert resolution_model.Q is None + assert len(resolution_model._component_collections) == 0 + + @pytest.mark.parametrize( + ('kwargs', 'match'), + [ + ( + {'sample_model': 'not_a_sample_model'}, + r'sample_model must be an instance of SampleModel', + ), + ( + {'normalize_area': 'not_a_bool'}, + r'normalize_area must be True or False', + ), + ( + {'fix_parameters': 'not_a_bool'}, + r'fix_parameters must be True or False', + ), + ], + ) + def test_from_sample_model_invalid_arguments( + self, + sample_model, + kwargs, + match, + ): + # WHEN + valid_kwargs = { + 'sample_model': sample_model, + 'normalize_area': False, + 'fix_parameters': False, + } + valid_kwargs.update(kwargs) + + # THEN EXPECT + with pytest.raises(TypeError, match=match): + ResolutionModel.from_sample_model( + **valid_kwargs, + ) + + def test_from_sample_model_invalid_components(self, sample_model): + # WHEN + invalid_component = DeltaFunction(name='InvalidDelta') + sample_model.append_component(invalid_component) + + # THEN EXPECT + with pytest.raises( + TypeError, + match='cannot be a DeltaFunction', + ): + ResolutionModel.from_sample_model(sample_model)