From 03f0da327b2d722e759cadd35bdd5a95dfb6c07e Mon Sep 17 00:00:00 2001 From: Piotr Rozyczko Date: Sat, 31 Jan 2026 11:52:22 +0100 Subject: [PATCH 1/9] bilayer class --- src/easyreflectometry/sample/__init__.py | 2 + .../sample/assemblies/bilayer.py | 425 ++++++++++++++++++ tests/sample/assemblies/test_bilayer.py | 410 +++++++++++++++++ 3 files changed, 837 insertions(+) create mode 100644 src/easyreflectometry/sample/assemblies/bilayer.py create mode 100644 tests/sample/assemblies/test_bilayer.py diff --git a/src/easyreflectometry/sample/__init__.py b/src/easyreflectometry/sample/__init__.py index 2012cd44..4991b975 100644 --- a/src/easyreflectometry/sample/__init__.py +++ b/src/easyreflectometry/sample/__init__.py @@ -1,4 +1,5 @@ from .assemblies.base_assembly import BaseAssembly +from .assemblies.bilayer import Bilayer from .assemblies.gradient_layer import GradientLayer from .assemblies.multilayer import Multilayer from .assemblies.repeating_multilayer import RepeatingMultilayer @@ -15,6 +16,7 @@ __all__ = ( 'BaseAssembly', + 'Bilayer', 'GradientLayer', 'Layer', 'LayerAreaPerMolecule', diff --git a/src/easyreflectometry/sample/assemblies/bilayer.py b/src/easyreflectometry/sample/assemblies/bilayer.py new file mode 100644 index 00000000..888c4d1d --- /dev/null +++ b/src/easyreflectometry/sample/assemblies/bilayer.py @@ -0,0 +1,425 @@ +from __future__ import annotations + +from typing import Optional + +from easyscience import global_object +from easyscience.variable import Parameter + +from ..collections.layer_collection import LayerCollection +from ..elements.layers.layer_area_per_molecule import LayerAreaPerMolecule +from ..elements.materials.material import Material +from .base_assembly import BaseAssembly + + +class Bilayer(BaseAssembly): + """A lipid bilayer consisting of two surfactant layers where one is inverted. + + The bilayer structure is: Front Head - Front Tail - Back Tail - Back Head + + This assembly comes pre-populated with physically meaningful constraints: + - Both tail layers share the same structural parameters (thickness, area per molecule) + - Head layers share thickness and area per molecule (different hydration/solvent_fraction allowed) + - A single roughness parameter applies to all interfaces (conformal roughness) + + More information about the usage of this assembly is available in the + `bilayer documentation`_ + + .. _`bilayer documentation`: ../sample/assemblies_library.html#bilayer + """ + + def __init__( + self, + front_head_layer: Optional[LayerAreaPerMolecule] = None, + tail_layer: Optional[LayerAreaPerMolecule] = None, + back_head_layer: Optional[LayerAreaPerMolecule] = None, + name: str = 'EasyBilayer', + unique_name: Optional[str] = None, + constrain_heads: bool = True, + conformal_roughness: bool = True, + interface=None, + ): + """Constructor. + + :param front_head_layer: Layer representing the front head part of the bilayer. + :param tail_layer: Layer representing the tail part of the bilayer. A second tail + layer is created internally with parameters constrained to this one. + :param back_head_layer: Layer representing the back head part of the bilayer. + :param name: Name for bilayer, defaults to 'EasyBilayer'. + :param constrain_heads: Constrain head layer thickness and area per molecule, defaults to `True`. + :param conformal_roughness: Constrain roughness to be the same for all layers, defaults to `True`. + :param interface: Calculator interface, defaults to `None`. + """ + # Generate unique name for nested objects + if unique_name is None: + unique_name = global_object.generate_unique_name(self.__class__.__name__) + + # Create default front head layer if not provided + if front_head_layer is None: + d2o_front = Material( + sld=6.36, + isld=0, + name='D2O', + unique_name=unique_name + '_MaterialFrontHead', + interface=interface, + ) + front_head_layer = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=10.0, + solvent=d2o_front, + solvent_fraction=0.2, + area_per_molecule=48.2, + roughness=3.0, + name='DPPC Head Front', + unique_name=unique_name + '_LayerAreaPerMoleculeFrontHead', + interface=interface, + ) + + # Create default tail layer if not provided + if tail_layer is None: + air = Material( + sld=0, + isld=0, + name='Air', + unique_name=unique_name + '_MaterialTail', + interface=interface, + ) + tail_layer = LayerAreaPerMolecule( + molecular_formula='C32D64', + thickness=16.0, + solvent=air, + solvent_fraction=0.0, + area_per_molecule=48.2, + roughness=3.0, + name='DPPC Tail', + unique_name=unique_name + '_LayerAreaPerMoleculeTail', + interface=interface, + ) + + # Create second tail layer with same parameters as the first + # This will be constrained to the first tail layer + air_back = Material( + sld=0, + isld=0, + name='Air', + unique_name=unique_name + '_MaterialBackTail', + interface=interface, + ) + back_tail_layer = LayerAreaPerMolecule( + molecular_formula=tail_layer.molecular_formula, + thickness=tail_layer.thickness.value, + solvent=air_back, + solvent_fraction=tail_layer.solvent_fraction, + area_per_molecule=tail_layer.area_per_molecule, + roughness=tail_layer.roughness.value, + name=tail_layer.name + ' Back', + unique_name=unique_name + '_LayerAreaPerMoleculeBackTail', + interface=interface, + ) + + # Create default back head layer if not provided + if back_head_layer is None: + d2o_back = Material( + sld=6.36, + isld=0, + name='D2O', + unique_name=unique_name + '_MaterialBackHead', + interface=interface, + ) + back_head_layer = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=10.0, + solvent=d2o_back, + solvent_fraction=0.2, + area_per_molecule=48.2, + roughness=3.0, + name='DPPC Head Back', + unique_name=unique_name + '_LayerAreaPerMoleculeBackHead', + interface=interface, + ) + + # Store the front tail layer reference for constraint setup + self._front_tail_layer = tail_layer + self._back_tail_layer = back_tail_layer + + # Create layer collection: front_head, front_tail, back_tail, back_head + bilayer_layers = LayerCollection( + front_head_layer, + tail_layer, + back_tail_layer, + back_head_layer, + name='Layers', + unique_name=unique_name + '_LayerCollection', + interface=interface, + ) + + super().__init__( + name=name, + unique_name=unique_name, + type='Bilayer', + layers=bilayer_layers, + interface=interface, + ) + + self.interface = interface + self._conformal_roughness = False + self._constrain_heads = False + self._tail_constraints_setup = False + + # Setup tail layer constraints (back tail depends on front tail) + self._setup_tail_constraints() + + # Apply head constraints if requested + if constrain_heads: + self.constrain_heads = True + + # Apply conformal roughness if requested + if conformal_roughness: + self.conformal_roughness = True + + def _setup_tail_constraints(self) -> None: + """Setup constraints so back tail layer parameters depend on front tail layer.""" + front_tail = self._front_tail_layer + back_tail = self._back_tail_layer + + # Constrain thickness + back_tail.thickness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_tail.thickness}, + ) + + # Constrain area per molecule + back_tail._area_per_molecule.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_tail._area_per_molecule}, + ) + + # Constrain solvent fraction + back_tail.material._fraction.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_tail.material._fraction}, + ) + + self._tail_constraints_setup = True + + @property + def front_head_layer(self) -> Optional[LayerAreaPerMolecule]: + """Get the front head layer of the bilayer.""" + return self.layers[0] + + @front_head_layer.setter + def front_head_layer(self, layer: LayerAreaPerMolecule) -> None: + """Set the front head layer of the bilayer.""" + self.layers[0] = layer + + @property + def front_tail_layer(self) -> Optional[LayerAreaPerMolecule]: + """Get the front tail layer of the bilayer.""" + return self.layers[1] + + @property + def tail_layer(self) -> Optional[LayerAreaPerMolecule]: + """Get the tail layer (alias for front_tail_layer for serialization compatibility).""" + return self.front_tail_layer + + @property + def back_tail_layer(self) -> Optional[LayerAreaPerMolecule]: + """Get the back tail layer of the bilayer.""" + return self.layers[2] + + @property + def back_head_layer(self) -> Optional[LayerAreaPerMolecule]: + """Get the back head layer of the bilayer.""" + return self.layers[3] + + @back_head_layer.setter + def back_head_layer(self, layer: LayerAreaPerMolecule) -> None: + """Set the back head layer of the bilayer.""" + self.layers[3] = layer + + @property + def constrain_heads(self) -> bool: + """Get the head layer constraint status.""" + return self._constrain_heads + + @constrain_heads.setter + def constrain_heads(self, status: bool) -> None: + """Set the status for head layer constraints. + + When enabled, the back head layer thickness and area per molecule + are constrained to match the front head layer. Solvent fraction + (hydration) remains independent. + + :param status: Boolean for the constraint status. + """ + if status: + self._enable_head_constraints() + else: + self._disable_head_constraints() + self._constrain_heads = status + + def _enable_head_constraints(self) -> None: + """Enable head layer constraints (thickness and area per molecule).""" + front_head = self.front_head_layer + back_head = self.back_head_layer + + # Constrain thickness + back_head.thickness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_head.thickness}, + ) + + # Constrain area per molecule + back_head._area_per_molecule.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_head._area_per_molecule}, + ) + + def _disable_head_constraints(self) -> None: + """Disable head layer constraints.""" + self.back_head_layer.thickness.make_independent() + self.back_head_layer._area_per_molecule.make_independent() + + @property + def conformal_roughness(self) -> bool: + """Get the roughness constraint status.""" + return self._conformal_roughness + + @conformal_roughness.setter + def conformal_roughness(self, status: bool) -> None: + """Set the status for conformal roughness. + + When enabled, all layers share the same roughness parameter + (controlled by the front head layer). + + :param status: Boolean for the constraint status. + """ + if status: + self._setup_roughness_constraints() + self._enable_roughness_constraints() + else: + if self._roughness_constraints_setup: + self._disable_roughness_constraints() + self._conformal_roughness = status + + def constrain_solvent_roughness(self, solvent_roughness: Parameter) -> None: + """Add the constraint to the solvent roughness. + + :param solvent_roughness: The solvent roughness parameter. + """ + if not self.conformal_roughness: + raise ValueError('Roughness must be conformal to use this function.') + solvent_roughness.value = self.front_head_layer.roughness.value + solvent_roughness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': self.front_head_layer.roughness}, + ) + + def constrain_multiple_contrast( + self, + another_contrast: Bilayer, + front_head_thickness: bool = True, + back_head_thickness: bool = True, + tail_thickness: bool = True, + front_head_area_per_molecule: bool = True, + back_head_area_per_molecule: bool = True, + tail_area_per_molecule: bool = True, + front_head_fraction: bool = True, + back_head_fraction: bool = True, + tail_fraction: bool = True, + ) -> None: + """Constrain structural parameters between bilayer objects. + + :param another_contrast: The bilayer to constrain to. + :param front_head_thickness: Constrain front head thickness. + :param back_head_thickness: Constrain back head thickness. + :param tail_thickness: Constrain tail thickness. + :param front_head_area_per_molecule: Constrain front head area per molecule. + :param back_head_area_per_molecule: Constrain back head area per molecule. + :param tail_area_per_molecule: Constrain tail area per molecule. + :param front_head_fraction: Constrain front head solvent fraction. + :param back_head_fraction: Constrain back head solvent fraction. + :param tail_fraction: Constrain tail solvent fraction. + """ + if front_head_thickness: + self.front_head_layer.thickness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_head_layer.thickness}, + ) + + if back_head_thickness: + self.back_head_layer.thickness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.back_head_layer.thickness}, + ) + + if tail_thickness: + self.front_tail_layer.thickness.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_tail_layer.thickness}, + ) + + if front_head_area_per_molecule: + self.front_head_layer._area_per_molecule.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_head_layer._area_per_molecule}, + ) + + if back_head_area_per_molecule: + self.back_head_layer._area_per_molecule.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.back_head_layer._area_per_molecule}, + ) + + if tail_area_per_molecule: + self.front_tail_layer._area_per_molecule.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_tail_layer._area_per_molecule}, + ) + + if front_head_fraction: + self.front_head_layer.material._fraction.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_head_layer.material._fraction}, + ) + + if back_head_fraction: + self.back_head_layer.material._fraction.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.back_head_layer.material._fraction}, + ) + + if tail_fraction: + self.front_tail_layer.material._fraction.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_tail_layer.material._fraction}, + ) + + @property + def _dict_repr(self) -> dict: + """A simplified dict representation.""" + return { + self.name: { + 'front_head_layer': self.front_head_layer._dict_repr, + 'front_tail_layer': self.front_tail_layer._dict_repr, + 'back_tail_layer': self.back_tail_layer._dict_repr, + 'back_head_layer': self.back_head_layer._dict_repr, + 'constrain_heads': self.constrain_heads, + 'conformal_roughness': self.conformal_roughness, + } + } + + def as_dict(self, skip: Optional[list[str]] = None) -> dict: + """Produce a cleaned dict using a custom as_dict method. + + The resulting dict matches the parameters in __init__ + + :param skip: List of keys to skip, defaults to `None`. + """ + this_dict = super().as_dict(skip=skip) + this_dict['front_head_layer'] = self.front_head_layer.as_dict(skip=skip) + this_dict['tail_layer'] = self.front_tail_layer.as_dict(skip=skip) + this_dict['back_head_layer'] = self.back_head_layer.as_dict(skip=skip) + this_dict['constrain_heads'] = self.constrain_heads + this_dict['conformal_roughness'] = self.conformal_roughness + del this_dict['layers'] + return this_dict diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py new file mode 100644 index 00000000..a17ae0b0 --- /dev/null +++ b/tests/sample/assemblies/test_bilayer.py @@ -0,0 +1,410 @@ +""" +Tests for Bilayer class module +""" + +__author__ = 'github.com/easyscience' +__version__ = '0.0.1' + + +from easyscience import global_object + +from easyreflectometry.sample.assemblies.bilayer import Bilayer +from easyreflectometry.sample.elements.layers.layer import Layer +from easyreflectometry.sample.elements.layers.layer_area_per_molecule import LayerAreaPerMolecule +from easyreflectometry.sample.elements.materials.material import Material + + +class TestBilayer: + def test_default(self): + """Test default bilayer creation with expected structure.""" + p = Bilayer() + assert p.name == 'EasyBilayer' + assert p._type == 'Bilayer' + + # Check layer count + assert len(p.layers) == 4 + + # Check layer order: front_head, front_tail, back_tail, back_head + assert p.layers[0].name == 'DPPC Head Front' + assert p.front_head_layer.name == 'DPPC Head Front' + + assert p.layers[1].name == 'DPPC Tail' + assert p.front_tail_layer.name == 'DPPC Tail' + + assert p.layers[2].name == 'DPPC Tail Back' + assert p.back_tail_layer.name == 'DPPC Tail Back' + + assert p.layers[3].name == 'DPPC Head Back' + assert p.back_head_layer.name == 'DPPC Head Back' + + def test_default_constraints_enabled(self): + """Test that default bilayer has constraints enabled.""" + p = Bilayer() + + # Default should have conformal roughness and head constraints + assert p.conformal_roughness is True + assert p.constrain_heads is True + + def test_layer_structure(self): + """Verify 4 layers in correct order.""" + p = Bilayer() + + assert p.front_head_layer is p.layers[0] + assert p.front_tail_layer is p.layers[1] + assert p.back_tail_layer is p.layers[2] + assert p.back_head_layer is p.layers[3] + + def test_custom_layers(self): + """Test creation with custom head/tail layers.""" + d2o = Material(sld=6.36, isld=0, name='D2O') + air = Material(sld=0, isld=0, name='Air') + + front_head = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=12.0, + solvent=d2o, + solvent_fraction=0.3, + area_per_molecule=50.0, + roughness=2.0, + name='Custom Front Head', + ) + tail = LayerAreaPerMolecule( + molecular_formula='C32D64', + thickness=18.0, + solvent=air, + solvent_fraction=0.0, + area_per_molecule=50.0, + roughness=2.0, + name='Custom Tail', + ) + back_head = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=12.0, + solvent=d2o, + solvent_fraction=0.4, # Different hydration + area_per_molecule=50.0, + roughness=2.0, + name='Custom Back Head', + ) + + p = Bilayer( + front_head_layer=front_head, + tail_layer=tail, + back_head_layer=back_head, + name='Custom Bilayer', + ) + + assert p.name == 'Custom Bilayer' + assert p.front_head_layer.name == 'Custom Front Head' + assert p.front_tail_layer.name == 'Custom Tail' + assert p.back_head_layer.name == 'Custom Back Head' + assert p.front_head_layer.thickness.value == 12.0 + assert p.front_tail_layer.thickness.value == 18.0 + + def test_tail_layers_linked(self): + """Test that both tail layers share parameters.""" + p = Bilayer() + + # Initial values should match + assert p.front_tail_layer.thickness.value == p.back_tail_layer.thickness.value + assert p.front_tail_layer.area_per_molecule == p.back_tail_layer.area_per_molecule + + # Change front tail thickness - back tail should follow + p.front_tail_layer.thickness.value = 20.0 + assert p.front_tail_layer.thickness.value == 20.0 + assert p.back_tail_layer.thickness.value == 20.0 + + # Change front tail area per molecule - back tail should follow + p.front_tail_layer._area_per_molecule.value = 55.0 + assert p.front_tail_layer.area_per_molecule == 55.0 + assert p.back_tail_layer.area_per_molecule == 55.0 + + def test_constrain_heads_enabled(self): + """Test head thickness/area constraint when enabled.""" + p = Bilayer(constrain_heads=True) + + # Change front head thickness - back head should follow + p.front_head_layer.thickness.value = 15.0 + assert p.front_head_layer.thickness.value == 15.0 + assert p.back_head_layer.thickness.value == 15.0 + + # Change front head area per molecule - back head should follow + p.front_head_layer._area_per_molecule.value = 60.0 + assert p.front_head_layer.area_per_molecule == 60.0 + assert p.back_head_layer.area_per_molecule == 60.0 + + def test_constrain_heads_disabled(self): + """Test heads are independent when constraint disabled.""" + p = Bilayer(constrain_heads=False) + + # Set different values for front and back heads + p.front_head_layer.thickness.value = 15.0 + p.back_head_layer.thickness.value = 12.0 + + assert p.front_head_layer.thickness.value == 15.0 + assert p.back_head_layer.thickness.value == 12.0 + + def test_constrain_heads_toggle(self): + """Test enabling/disabling head constraints at runtime.""" + p = Bilayer(constrain_heads=False) + + # Set different values + p.front_head_layer.thickness.value = 15.0 + p.back_head_layer.thickness.value = 12.0 + + # Enable constraint - back head should match front head + p.constrain_heads = True + assert p.constrain_heads is True + + # Change front head - back should follow + p.front_head_layer.thickness.value = 20.0 + assert p.back_head_layer.thickness.value == 20.0 + + # Disable constraint + p.constrain_heads = False + assert p.constrain_heads is False + + # Now they can be independent + p.back_head_layer.thickness.value = 18.0 + assert p.front_head_layer.thickness.value == 20.0 + assert p.back_head_layer.thickness.value == 18.0 + + def test_head_hydration_independent(self): + """Test that head hydrations remain independent even with constraints.""" + p = Bilayer(constrain_heads=True) + + # Set different solvent fractions + p.front_head_layer.solvent_fraction = 0.3 + p.back_head_layer.solvent_fraction = 0.5 + + # They should remain independent + assert p.front_head_layer.solvent_fraction == 0.3 + assert p.back_head_layer.solvent_fraction == 0.5 + + def test_conformal_roughness_enabled(self): + """Test all roughnesses are linked when conformal roughness enabled.""" + p = Bilayer(conformal_roughness=True) + + # Change front head roughness - all should follow + p.front_head_layer.roughness.value = 5.0 + assert p.front_head_layer.roughness.value == 5.0 + assert p.front_tail_layer.roughness.value == 5.0 + assert p.back_tail_layer.roughness.value == 5.0 + assert p.back_head_layer.roughness.value == 5.0 + + def test_conformal_roughness_disabled(self): + """Test roughnesses are independent when conformal roughness disabled.""" + p = Bilayer(conformal_roughness=False) + + # Set different roughnesses + p.front_head_layer.roughness.value = 2.0 + p.front_tail_layer.roughness.value = 3.0 + p.back_tail_layer.roughness.value = 4.0 + p.back_head_layer.roughness.value = 5.0 + + assert p.front_head_layer.roughness.value == 2.0 + assert p.front_tail_layer.roughness.value == 3.0 + assert p.back_tail_layer.roughness.value == 4.0 + assert p.back_head_layer.roughness.value == 5.0 + + def test_conformal_roughness_toggle(self): + """Test enabling/disabling conformal roughness at runtime.""" + p = Bilayer(conformal_roughness=False) + + # Set different values + p.front_head_layer.roughness.value = 2.0 + p.back_head_layer.roughness.value = 5.0 + + # Enable conformal roughness + p.conformal_roughness = True + assert p.conformal_roughness is True + + # Change front head - all should follow + p.front_head_layer.roughness.value = 4.0 + assert p.front_tail_layer.roughness.value == 4.0 + assert p.back_tail_layer.roughness.value == 4.0 + assert p.back_head_layer.roughness.value == 4.0 + + # Disable conformal roughness + p.conformal_roughness = False + assert p.conformal_roughness is False + + def test_get_set_front_head_layer(self): + """Test getting and setting front head layer.""" + p = Bilayer() + new_layer = LayerAreaPerMolecule( + molecular_formula='C8H16NO6P', + thickness=8.0, + name='New Front Head', + ) + + p.front_head_layer = new_layer + + assert p.front_head_layer == new_layer + assert p.layers[0] == new_layer + + def test_get_set_back_head_layer(self): + """Test getting and setting back head layer.""" + p = Bilayer() + new_layer = LayerAreaPerMolecule( + molecular_formula='C8H16NO6P', + thickness=8.0, + name='New Back Head', + ) + + p.back_head_layer = new_layer + + assert p.back_head_layer == new_layer + assert p.layers[3] == new_layer + + def test_dict_repr(self): + """Test dictionary representation.""" + p = Bilayer() + + dict_repr = p._dict_repr + assert 'EasyBilayer' in dict_repr + assert 'front_head_layer' in dict_repr['EasyBilayer'] + assert 'front_tail_layer' in dict_repr['EasyBilayer'] + assert 'back_tail_layer' in dict_repr['EasyBilayer'] + assert 'back_head_layer' in dict_repr['EasyBilayer'] + assert 'constrain_heads' in dict_repr['EasyBilayer'] + assert 'conformal_roughness' in dict_repr['EasyBilayer'] + + +def test_dict_round_trip(): + """Test serialization/deserialization round trip.""" + # When + d2o = Material(sld=6.36, isld=0, name='D2O') + air = Material(sld=0, isld=0, name='Air') + + front_head = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=12.0, + solvent=d2o, + solvent_fraction=0.3, + area_per_molecule=50.0, + roughness=2.0, + name='Custom Front Head', + ) + tail = LayerAreaPerMolecule( + molecular_formula='C32D64', + thickness=18.0, + solvent=air, + solvent_fraction=0.0, + area_per_molecule=50.0, + roughness=2.0, + name='Custom Tail', + ) + back_head = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=12.0, + solvent=d2o, + solvent_fraction=0.4, + area_per_molecule=50.0, + roughness=2.0, + name='Custom Back Head', + ) + + p = Bilayer( + front_head_layer=front_head, + tail_layer=tail, + back_head_layer=back_head, + constrain_heads=False, + conformal_roughness=False, + ) + p_dict = p.as_dict() + global_object.map._clear() + + # Then + q = Bilayer.from_dict(p_dict) + + # Expect + assert sorted(p.as_dict()) == sorted(q.as_dict()) + + +def test_dict_round_trip_constraints_enabled(): + """Test round trip with constraints enabled.""" + # When + p = Bilayer(constrain_heads=True, conformal_roughness=True) + p_dict = p.as_dict() + global_object.map._clear() + + # Then + q = Bilayer.from_dict(p_dict) + + # Expect + assert q.constrain_heads is True + assert q.conformal_roughness is True + assert sorted(p.as_dict()) == sorted(q.as_dict()) + + +def test_dict_round_trip_constraints_disabled(): + """Test round trip with constraints disabled.""" + # When + p = Bilayer(constrain_heads=False, conformal_roughness=False) + p_dict = p.as_dict() + global_object.map._clear() + + # Then + q = Bilayer.from_dict(p_dict) + + # Expect + assert q.constrain_heads is False + assert q.conformal_roughness is False + assert sorted(p.as_dict()) == sorted(q.as_dict()) + + +def test_constrain_multiple_contrast(): + """Test multi-contrast constraints between bilayers.""" + # When + p1 = Bilayer(name='Bilayer 1', constrain_heads=False) + p2 = Bilayer(name='Bilayer 2', constrain_heads=False) + + # Set initial values + p1.front_head_layer.thickness.value = 10.0 + p1.front_tail_layer.thickness.value = 16.0 + + # Constrain p2 to p1 + p2.constrain_multiple_contrast( + p1, + front_head_thickness=True, + tail_thickness=True, + ) + + # Then - p2 values should match p1 + assert p2.front_head_layer.thickness.value == 10.0 + assert p2.front_tail_layer.thickness.value == 16.0 + + # Change p1 - p2 should follow + p1.front_head_layer.thickness.value = 12.0 + assert p2.front_head_layer.thickness.value == 12.0 + + +def test_constrain_solvent_roughness(): + """Test constraining solvent roughness to bilayer roughness.""" + # When + p = Bilayer(conformal_roughness=True) + layer = Layer() + + p.front_head_layer.roughness.value = 4.0 + + # Then + p.constrain_solvent_roughness(layer.roughness) + + # Expect + assert layer.roughness.value == 4.0 + + # Change bilayer roughness - solvent should follow + p.front_head_layer.roughness.value = 5.0 + assert layer.roughness.value == 5.0 + + +def test_constrain_solvent_roughness_error(): + """Test error when constraining solvent roughness without conformal roughness.""" + import pytest + + p = Bilayer(conformal_roughness=False) + layer = Layer() + + with pytest.raises(ValueError, match='Roughness must be conformal'): + p.constrain_solvent_roughness(layer.roughness) From fac389ff4f6baee86ccb91e1aeb2caad62dbeb80 Mon Sep 17 00:00:00 2001 From: Piotr Rozyczko Date: Sat, 31 Jan 2026 14:52:57 +0100 Subject: [PATCH 2/9] bilayer notebook and docs --- docs/src/api/project.rst | 7 + .../tutorials/basic/assemblies_library.rst | 91 ++- docs/src/tutorials/simulation/bilayer.ipynb | 542 ++++++++++++++++++ docs/src/tutorials/simulation/simulation.rst | 1 + 4 files changed, 640 insertions(+), 1 deletion(-) create mode 100644 docs/src/api/project.rst create mode 100644 docs/src/tutorials/simulation/bilayer.ipynb diff --git a/docs/src/api/project.rst b/docs/src/api/project.rst new file mode 100644 index 00000000..277157b0 --- /dev/null +++ b/docs/src/api/project.rst @@ -0,0 +1,7 @@ +Project +======= + +.. autoclass:: easyreflectometry.project.Project + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/src/tutorials/basic/assemblies_library.rst b/docs/src/tutorials/basic/assemblies_library.rst index 7e79844f..d72201a9 100644 --- a/docs/src/tutorials/basic/assemblies_library.rst +++ b/docs/src/tutorials/basic/assemblies_library.rst @@ -153,8 +153,97 @@ Furthermore, as shown in the `surfactant monolayer tutorial`_ the conformal roug The use of the :py:class:`SurfactantLayer` in multiple contrast data analysis is shown in a `multiple contrast tutorial`_. +:py:class:`Bilayer` +------------------- + +The :py:class:`Bilayer` assembly type represents a phospholipid bilayer at an interface. +It consists of two surfactant layers where one is inverted, creating the structure: + +.. code-block:: text + + Head₁ - Tail₁ - Tail₂ - Head₂ + +This assembly is particularly useful for studying supported lipid bilayers and membrane systems. +The bilayer comes pre-populated with physically meaningful constraints: + +- Both tail layers share the same structural parameters (thickness, area per molecule) +- Head layers share thickness and area per molecule (different hydration/solvent fraction allowed) +- A single roughness parameter applies to all interfaces (conformal roughness) + +These default constraints can be enabled or disabled as needed for specific analyses. + +The creation of a :py:class:`Bilayer` object is shown below. + +.. code-block:: python + + from easyreflectometry.sample import Bilayer + from easyreflectometry.sample import LayerAreaPerMolecule + from easyreflectometry.sample import Material + + # Create materials for solvents + d2o = Material(sld=6.36, isld=0.0, name='D2O') + air = Material(sld=0.0, isld=0.0, name='Air') + + # Create head layer (used for front, back head will be auto-created with constraints) + head = LayerAreaPerMolecule( + molecular_formula='C10H18NO8P', + thickness=10.0, + solvent=d2o, + solvent_fraction=0.3, + area_per_molecule=48.2, + roughness=3.0, + name='DPPC Head' + ) + + # Create tail layer (both tail positions will share these parameters) + tail = LayerAreaPerMolecule( + molecular_formula='C32D64', + thickness=16.0, + solvent=air, + solvent_fraction=0.0, + area_per_molecule=48.2, + roughness=3.0, + name='DPPC Tail' + ) + + # Create bilayer with default constraints + bilayer = Bilayer( + front_head_layer=head, + tail_layer=tail, + constrain_heads=True, + conformal_roughness=True, + name='DPPC Bilayer' + ) + +The head layers can have different solvent fractions (hydration) even when constrained, +enabling the modeling of asymmetric bilayers at interfaces where the two sides of the +bilayer may have different solvent exposure. + +The constraints can be controlled at runtime: + +.. code-block:: python + + # Disable head constraints to allow different head layer structures + bilayer.constrain_heads = False + + # Disable conformal roughness to allow different roughness values + bilayer.conformal_roughness = False + +Individual layers can be accessed via properties: + +.. code-block:: python + + # Access the four layers + bilayer.front_head_layer # First head layer + bilayer.front_tail_layer # First tail layer + bilayer.back_tail_layer # Second tail layer (constrained to front tail) + bilayer.back_head_layer # Second head layer + +For more detailed examples including simulation and parameter access, see the `bilayer tutorial`_. + .. _`simple fitting tutorial`: ../tutorials/simple_fitting.html .. _`tutorial`: ../tutorials/repeating.html .. _`surfactant monolayer tutorial`: ../tutorials/monolayer.html -.. _`multiple contrast tutorial`: ../tutorials/multi_contrast.html \ No newline at end of file +.. _`multiple contrast tutorial`: ../tutorials/multi_contrast.html +.. _`bilayer tutorial`: ../tutorials/simulation/bilayer.html \ No newline at end of file diff --git a/docs/src/tutorials/simulation/bilayer.ipynb b/docs/src/tutorials/simulation/bilayer.ipynb new file mode 100644 index 00000000..9725ab1e --- /dev/null +++ b/docs/src/tutorials/simulation/bilayer.ipynb @@ -0,0 +1,542 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b5afc8c4", + "metadata": {}, + "source": [ + "# Simulating a Phospholipid Bilayer\n", + "\n", + "Phospholipid bilayers are fundamental structures in biological membranes and are commonly studied using neutron and X-ray reflectometry.\n", + "In this tutorial, we will explore how to use the `Bilayer` assembly in `easyreflectometry` to model a lipid bilayer structure.\n", + "\n", + "A bilayer consists of two surfactant layers arranged in an inverted configuration:\n", + "\n", + "```\n", + "Head₁ - Tail₁ - Tail₂ - Head₂\n", + "```\n", + "\n", + "The `Bilayer` assembly comes with pre-configured constraints that represent physically meaningful relationships:\n", + "- Both tail layers share the same structural parameters\n", + "- Head layers share thickness and area per molecule (but can have different hydration)\n", + "- Conformal roughness across all interfaces" + ] + }, + { + "cell_type": "markdown", + "id": "f6021005", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "First, we import the necessary modules and configure matplotlib for inline plotting." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1bfa00f4", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import easyreflectometry\n", + "from easyreflectometry.calculators import CalculatorFactory\n", + "from easyreflectometry.sample import Bilayer\n", + "from easyreflectometry.sample import LayerAreaPerMolecule\n", + "from easyreflectometry.sample import Material\n", + "from easyreflectometry.sample import Layer\n", + "from easyreflectometry.sample import Multilayer\n", + "from easyreflectometry.sample import Sample\n", + "from easyreflectometry.model import Model\n", + "from easyreflectometry.model import PercentageFwhm" + ] + }, + { + "cell_type": "markdown", + "id": "be41f6f6", + "metadata": {}, + "source": [ + "## Creating a Bilayer\n", + "\n", + "We'll create a DPPC (dipalmitoylphosphatidylcholine) bilayer, a common model phospholipid.\n", + "\n", + "First, let's define the materials for our solvents - D₂O (heavy water) for the head groups and air for the tail groups." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76bd9056", + "metadata": {}, + "outputs": [], + "source": [ + "# Define solvent materials\n", + "d2o = Material(sld=6.36, isld=0, name='D2O')\n", + "air = Material(sld=0, isld=0, name='Air')" + ] + }, + { + "cell_type": "markdown", + "id": "34a0da6c", + "metadata": {}, + "source": [ + "### Creating Layer Components\n", + "\n", + "Now we create the head and tail layers using `LayerAreaPerMolecule`. This approach allows us to define layers based on their chemical formula and area per molecule, which provides a more physically meaningful parameterization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54980c8e", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a head layer for the bilayer\n", + "# The head group of DPPC has formula C10H18NO8P\n", + "head_layer = LayerAreaPerMolecule(\n", + " molecular_formula='C10H18NO8P',\n", + " thickness=10.0,\n", + " solvent=d2o,\n", + " solvent_fraction=0.3, # 30% solvent in head region\n", + " area_per_molecule=48.2,\n", + " roughness=3.0,\n", + " name='DPPC Head'\n", + ")\n", + "\n", + "# Create a tail layer for the bilayer\n", + "# The tail group of deuterated DPPC has formula C32D64\n", + "tail_layer = LayerAreaPerMolecule(\n", + " molecular_formula='C32D64',\n", + " thickness=16.0,\n", + " solvent=air,\n", + " solvent_fraction=0.0, # No solvent in the tail region\n", + " area_per_molecule=48.2,\n", + " roughness=3.0,\n", + " name='DPPC Tail'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ce3083b9", + "metadata": {}, + "source": [ + "### Creating the Bilayer Assembly\n", + "\n", + "Now we create the `Bilayer` assembly. The bilayer will automatically:\n", + "- Create a second tail layer with parameters constrained to the first\n", + "- Create a back head layer with thickness and area per molecule constrained to the front head\n", + "- Apply conformal roughness across all layers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4afdd40f", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the bilayer with default constraints\n", + "bilayer = Bilayer(\n", + " front_head_layer=head_layer,\n", + " tail_layer=tail_layer,\n", + " constrain_heads=True, # Head layers share thickness and area per molecule\n", + " conformal_roughness=True, # All layers share the same roughness\n", + " name='DPPC Bilayer'\n", + ")\n", + "\n", + "print(bilayer)" + ] + }, + { + "cell_type": "markdown", + "id": "87e39d39", + "metadata": {}, + "source": [ + "## Exploring the Bilayer Structure\n", + "\n", + "Let's examine the layer structure of our bilayer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc5ffe80", + "metadata": {}, + "outputs": [], + "source": [ + "# The bilayer has 4 layers\n", + "print(f'Number of layers: {len(bilayer.layers)}')\n", + "print()\n", + "\n", + "# Access individual layers\n", + "print('Layer structure (from front to back):')\n", + "print(f' 1. Front Head: {bilayer.front_head_layer.name}')\n", + "print(f' 2. Front Tail: {bilayer.front_tail_layer.name}')\n", + "print(f' 3. Back Tail: {bilayer.back_tail_layer.name}')\n", + "print(f' 4. Back Head: {bilayer.back_head_layer.name}')" + ] + }, + { + "cell_type": "markdown", + "id": "d9b60931", + "metadata": {}, + "source": [ + "### Verifying Constraints\n", + "\n", + "Let's verify that the constraints are working correctly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a607ca1d", + "metadata": {}, + "outputs": [], + "source": [ + "# Check tail layer constraints - both tails should have the same parameters\n", + "print('Tail layer parameters:')\n", + "print(f' Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", + "print(f' Back tail thickness: {bilayer.back_tail_layer.thickness.value:.2f} Å')\n", + "print(f' Front tail APM: {bilayer.front_tail_layer.area_per_molecule:.2f} Ų')\n", + "print(f' Back tail APM: {bilayer.back_tail_layer.area_per_molecule:.2f} Ų')\n", + "print()\n", + "\n", + "# Change front tail thickness - back should follow\n", + "bilayer.front_tail_layer.thickness.value = 18.0\n", + "print('After changing front tail thickness to 18 Å:')\n", + "print(f' Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", + "print(f' Back tail thickness: {bilayer.back_tail_layer.thickness.value:.2f} Å')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffb84f80", + "metadata": {}, + "outputs": [], + "source": [ + "# Check head layer constraints\n", + "print('Head layer parameters:')\n", + "print(f' Front head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n", + "print(f' Back head thickness: {bilayer.back_head_layer.thickness.value:.2f} Å')\n", + "print()\n", + "\n", + "# Hydration (solvent fraction) remains independent\n", + "print('Head layer hydration (independent):')\n", + "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", + "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')\n", + "\n", + "# Change back head hydration - front should NOT change\n", + "bilayer.back_head_layer.solvent_fraction = 0.5\n", + "print()\n", + "print('After changing back head solvent fraction to 0.5:')\n", + "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", + "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39802839", + "metadata": {}, + "outputs": [], + "source": [ + "# Check conformal roughness\n", + "print('Roughness values (conformal):')\n", + "print(f' Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", + "print(f' Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back tail roughness: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back head roughness: {bilayer.back_head_layer.roughness.value:.2f} Å')\n", + "print()\n", + "\n", + "# Change front head roughness - all should follow\n", + "bilayer.front_head_layer.roughness.value = 4.0\n", + "print('After changing front head roughness to 4.0 Å:')\n", + "print(f' Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", + "print(f' Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back tail roughness: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back head roughness: {bilayer.back_head_layer.roughness.value:.2f} Å')" + ] + }, + { + "cell_type": "markdown", + "id": "84f1206b", + "metadata": {}, + "source": [ + "## Building a Complete Sample\n", + "\n", + "To simulate reflectometry, we need to create a complete sample with sub- and super-phases." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e509350", + "metadata": {}, + "outputs": [], + "source": [ + "# Reset bilayer parameters\n", + "bilayer.front_head_layer.roughness.value = 3.0\n", + "bilayer.front_tail_layer.thickness.value = 16.0\n", + "bilayer.back_head_layer.solvent_fraction = 0.3\n", + "\n", + "# Create superphase (air) and subphase (silicon substrate with oxide layer)\n", + "vacuum = Material(sld=0, isld=0, name='Vacuum')\n", + "superphase = Layer(material=vacuum, thickness=0, roughness=0, name='Vacuum Superphase')\n", + "\n", + "si = Material(sld=2.047, isld=0, name='Si')\n", + "sio2 = Material(sld=3.47, isld=0, name='SiO2')\n", + "si_layer = Layer(material=si, thickness=0, roughness=0, name='Si')\n", + "sio2_layer = Layer(material=sio2, thickness=15, roughness=3, name='SiO2')\n", + "\n", + "# D2O subphase\n", + "d2o_layer = Layer(material=d2o, thickness=0, roughness=3, name='D2O Subphase')\n", + "\n", + "# Create sample: superphase -> bilayer -> SiO2 -> Si\n", + "sample = Sample(\n", + " Multilayer(superphase, name='Superphase'),\n", + " bilayer,\n", + " Multilayer(d2o_layer, name='D2O'),\n", + " Multilayer([sio2_layer, si_layer], name='Substrate'),\n", + " name='Bilayer on Si/SiO2'\n", + ")\n", + "\n", + "print(sample)" + ] + }, + { + "cell_type": "markdown", + "id": "0ce27fe5", + "metadata": {}, + "source": [ + "## Simulating Reflectivity\n", + "\n", + "Now we can simulate the reflectometry profile for our bilayer sample." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58538a77", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the model\n", + "model = Model(\n", + " sample=sample,\n", + " scale=1.0,\n", + " background=1e-7,\n", + " resolution_function=PercentageFwhm(5),\n", + " name='Bilayer Model'\n", + ")\n", + "\n", + "# Set up the calculator\n", + "interface = CalculatorFactory()\n", + "model.interface = interface\n", + "\n", + "# Generate Q values\n", + "q = np.linspace(0.005, 0.3, 500)\n", + "\n", + "# Calculate reflectometry\n", + "reflectivity = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(10, 6))\n", + "plt.semilogy(q, reflectivity, 'b-', linewidth=2, label='Bilayer')\n", + "plt.xlabel('Q (Å⁻¹)')\n", + "plt.ylabel('Reflectivity')\n", + "plt.title('Simulated Reflectometry of DPPC Bilayer')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "06723f8f", + "metadata": {}, + "source": [ + "## Scattering Length Density Profile\n", + "\n", + "Let's also visualize the SLD profile of our bilayer structure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a61981e0", + "metadata": {}, + "outputs": [], + "source": [ + "# Get SLD profile\n", + "z, sld = model.interface().sld_profile(model.unique_name)\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.plot(z, sld, 'b-', linewidth=2)\n", + "plt.xlabel('Distance from interface (Å)')\n", + "plt.ylabel('SLD (10⁻⁶ Å⁻²)')\n", + "plt.title('SLD Profile of DPPC Bilayer on Si/SiO2')\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "504f5fbe", + "metadata": {}, + "source": [ + "## Modifying Constraints\n", + "\n", + "The bilayer constraints can be modified at runtime. Let's see how disabling conformal roughness affects the model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12bffbbb", + "metadata": {}, + "outputs": [], + "source": [ + "# Disable conformal roughness\n", + "bilayer.conformal_roughness = False\n", + "\n", + "# Now we can set different roughness values for each layer\n", + "bilayer.front_head_layer.roughness.value = 2.0\n", + "bilayer.front_tail_layer.roughness.value = 1.5\n", + "bilayer.back_tail_layer.roughness.value = 1.5\n", + "bilayer.back_head_layer.roughness.value = 4.0\n", + "\n", + "print('Individual roughness values after disabling conformal roughness:')\n", + "print(f' Front head: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", + "print(f' Front tail: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back tail: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", + "print(f' Back head: {bilayer.back_head_layer.roughness.value:.2f} Å')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5fd26d95", + "metadata": {}, + "outputs": [], + "source": [ + "# Compare reflectometry with different roughness configurations\n", + "reflectivity_variable_roughness = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "# Reset to conformal roughness\n", + "bilayer.conformal_roughness = True\n", + "bilayer.front_head_layer.roughness.value = 3.0\n", + "reflectivity_conformal = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.semilogy(q, reflectivity_conformal, 'b-', linewidth=2, label='Conformal roughness (3.0 Å)')\n", + "plt.semilogy(q, reflectivity_variable_roughness, 'r--', linewidth=2, label='Variable roughness')\n", + "plt.xlabel('Q (Å⁻¹)')\n", + "plt.ylabel('Reflectivity')\n", + "plt.title('Effect of Roughness Configuration on Reflectometry')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "955e1320", + "metadata": {}, + "source": [ + "## Asymmetric Hydration\n", + "\n", + "One of the key features of the `Bilayer` assembly is the ability to model asymmetric hydration - where the two sides of the bilayer have different solvent exposure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "993ad6d4", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a bilayer with asymmetric hydration\n", + "# This is common in supported bilayers where one side is in contact with a substrate\n", + "\n", + "# Low hydration on front (substrate side)\n", + "bilayer.front_head_layer.solvent_fraction = 0.1\n", + "\n", + "# Higher hydration on back (solution side)\n", + "bilayer.back_head_layer.solvent_fraction = 0.4\n", + "\n", + "print('Asymmetric hydration:')\n", + "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", + "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')\n", + "\n", + "# Calculate new reflectometry\n", + "reflectivity_asymmetric = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "# Reset to symmetric for comparison\n", + "bilayer.front_head_layer.solvent_fraction = 0.3\n", + "bilayer.back_head_layer.solvent_fraction = 0.3\n", + "reflectivity_symmetric = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.semilogy(q, reflectivity_symmetric, 'b-', linewidth=2, label='Symmetric hydration (0.3/0.3)')\n", + "plt.semilogy(q, reflectivity_asymmetric, 'r--', linewidth=2, label='Asymmetric hydration (0.1/0.4)')\n", + "plt.xlabel('Q (Å⁻¹)')\n", + "plt.ylabel('Reflectivity')\n", + "plt.title('Effect of Asymmetric Hydration on Reflectometry')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "df3d351f", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "In this tutorial, we have explored the `Bilayer` assembly in `easyreflectometry`:\n", + "\n", + "1. **Creating a bilayer**: Using `LayerAreaPerMolecule` components for head and tail layers\n", + "2. **Understanding constraints**: \n", + " - Tail layers share all structural parameters\n", + " - Head layers share thickness and area per molecule (hydration is independent)\n", + " - Conformal roughness applies to all layers by default\n", + "3. **Building a sample**: Combining the bilayer with sub- and super-phases\n", + "4. **Simulating reflectometry**: Using the calculator interface to generate reflectivity profiles\n", + "5. **Modifying constraints**: Disabling conformal roughness for more complex models\n", + "6. **Asymmetric hydration**: Modeling supported bilayers with different solvent exposure on each side\n", + "\n", + "The `Bilayer` assembly provides a convenient way to model phospholipid bilayers with physically meaningful constraints, reducing the number of free parameters while maintaining flexibility for complex systems." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "era", + "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.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/src/tutorials/simulation/simulation.rst b/docs/src/tutorials/simulation/simulation.rst index 4f187b05..5cf9ced4 100644 --- a/docs/src/tutorials/simulation/simulation.rst +++ b/docs/src/tutorials/simulation/simulation.rst @@ -6,5 +6,6 @@ These are basic simulation examples using the :py:mod:`easyreflectometry` librar .. toctree:: :maxdepth: 1 + bilayer.ipynb magnetism.ipynb resolution_functions.ipynb \ No newline at end of file From efe523ac21e148ec550a2058946c5be5eee3404b Mon Sep 17 00:00:00 2001 From: rozyczko Date: Mon, 2 Feb 2026 14:14:00 +0100 Subject: [PATCH 3/9] PR review #1 --- docs/src/tutorials/simulation/bilayer.ipynb | 10 ++-- pyproject.toml | 2 +- .../sample/assemblies/bilayer.py | 49 +++++++++---------- tests/sample/assemblies/test_bilayer.py | 12 ++--- 4 files changed, 34 insertions(+), 39 deletions(-) diff --git a/docs/src/tutorials/simulation/bilayer.ipynb b/docs/src/tutorials/simulation/bilayer.ipynb index 9725ab1e..5e8d3a0f 100644 --- a/docs/src/tutorials/simulation/bilayer.ipynb +++ b/docs/src/tutorials/simulation/bilayer.ipynb @@ -77,7 +77,7 @@ "source": [ "# Define solvent materials\n", "d2o = Material(sld=6.36, isld=0, name='D2O')\n", - "air = Material(sld=0, isld=0, name='Air')" + "air_matched_water = Material(sld=0, isld=0, name='Air Matched Water')" ] }, { @@ -111,10 +111,10 @@ "\n", "# Create a tail layer for the bilayer\n", "# The tail group of deuterated DPPC has formula C32D64\n", - "tail_layer = LayerAreaPerMolecule(\n", + "front_tail_layer = LayerAreaPerMolecule(\n", " molecular_formula='C32D64',\n", " thickness=16.0,\n", - " solvent=air,\n", + " solvent=air_matched_water,\n", " solvent_fraction=0.0, # No solvent in the tail region\n", " area_per_molecule=48.2,\n", " roughness=3.0,\n", @@ -145,7 +145,7 @@ "# Create the bilayer with default constraints\n", "bilayer = Bilayer(\n", " front_head_layer=head_layer,\n", - " tail_layer=tail_layer,\n", + " front_front_front_tail_layer=front_tail_layer,\n", " constrain_heads=True, # Head layers share thickness and area per molecule\n", " conformal_roughness=True, # All layers share the same roughness\n", " name='DPPC Bilayer'\n", @@ -534,7 +534,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.11" + "version": "3.12.12" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 630422ca..b618dc16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,7 @@ dependencies = [ "refnx", "refl1d>=1.0.0", "orsopy", - "svglib<1.6 ; platform_system=='Linux'", + "svglib<1.6 ; platform_system=='Linux' or sys_platform == 'darwin'", "xhtml2pdf", "bumps", ] diff --git a/src/easyreflectometry/sample/assemblies/bilayer.py b/src/easyreflectometry/sample/assemblies/bilayer.py index 888c4d1d..797ca792 100644 --- a/src/easyreflectometry/sample/assemblies/bilayer.py +++ b/src/easyreflectometry/sample/assemblies/bilayer.py @@ -30,7 +30,7 @@ class Bilayer(BaseAssembly): def __init__( self, front_head_layer: Optional[LayerAreaPerMolecule] = None, - tail_layer: Optional[LayerAreaPerMolecule] = None, + front_tail_layer: Optional[LayerAreaPerMolecule] = None, back_head_layer: Optional[LayerAreaPerMolecule] = None, name: str = 'EasyBilayer', unique_name: Optional[str] = None, @@ -41,7 +41,7 @@ def __init__( """Constructor. :param front_head_layer: Layer representing the front head part of the bilayer. - :param tail_layer: Layer representing the tail part of the bilayer. A second tail + :param front_tail_layer: Layer representing the tail part of the bilayer. A second tail layer is created internally with parameters constrained to this one. :param back_head_layer: Layer representing the back head part of the bilayer. :param name: Name for bilayer, defaults to 'EasyBilayer'. @@ -74,19 +74,19 @@ def __init__( interface=interface, ) - # Create default tail layer if not provided - if tail_layer is None: - air = Material( - sld=0, + # Create default front tail layer if not provided + if front_tail_layer is None: + d2o_tail = Material( + sld=6.36, isld=0, - name='Air', + name='D2O', unique_name=unique_name + '_MaterialTail', interface=interface, ) - tail_layer = LayerAreaPerMolecule( + front_tail_layer = LayerAreaPerMolecule( molecular_formula='C32D64', thickness=16.0, - solvent=air, + solvent=d2o_tail, solvent_fraction=0.0, area_per_molecule=48.2, roughness=3.0, @@ -97,21 +97,21 @@ def __init__( # Create second tail layer with same parameters as the first # This will be constrained to the first tail layer - air_back = Material( - sld=0, + d2o_back_tail = Material( + sld=6.36, isld=0, - name='Air', + name='D2O', unique_name=unique_name + '_MaterialBackTail', interface=interface, ) back_tail_layer = LayerAreaPerMolecule( - molecular_formula=tail_layer.molecular_formula, - thickness=tail_layer.thickness.value, - solvent=air_back, - solvent_fraction=tail_layer.solvent_fraction, - area_per_molecule=tail_layer.area_per_molecule, - roughness=tail_layer.roughness.value, - name=tail_layer.name + ' Back', + molecular_formula=front_tail_layer.molecular_formula, + thickness=front_tail_layer.thickness.value, + solvent=d2o_back_tail, + solvent_fraction=front_tail_layer.solvent_fraction, + area_per_molecule=front_tail_layer.area_per_molecule, + roughness=front_tail_layer.roughness.value, + name=front_tail_layer.name + ' Back', unique_name=unique_name + '_LayerAreaPerMoleculeBackTail', interface=interface, ) @@ -138,13 +138,13 @@ def __init__( ) # Store the front tail layer reference for constraint setup - self._front_tail_layer = tail_layer + self._front_tail_layer = front_tail_layer self._back_tail_layer = back_tail_layer # Create layer collection: front_head, front_tail, back_tail, back_head bilayer_layers = LayerCollection( front_head_layer, - tail_layer, + front_tail_layer, back_tail_layer, back_head_layer, name='Layers', @@ -216,11 +216,6 @@ def front_tail_layer(self) -> Optional[LayerAreaPerMolecule]: """Get the front tail layer of the bilayer.""" return self.layers[1] - @property - def tail_layer(self) -> Optional[LayerAreaPerMolecule]: - """Get the tail layer (alias for front_tail_layer for serialization compatibility).""" - return self.front_tail_layer - @property def back_tail_layer(self) -> Optional[LayerAreaPerMolecule]: """Get the back tail layer of the bilayer.""" @@ -417,7 +412,7 @@ def as_dict(self, skip: Optional[list[str]] = None) -> dict: """ this_dict = super().as_dict(skip=skip) this_dict['front_head_layer'] = self.front_head_layer.as_dict(skip=skip) - this_dict['tail_layer'] = self.front_tail_layer.as_dict(skip=skip) + this_dict['front_tail_layer'] = self.front_tail_layer.as_dict(skip=skip) this_dict['back_head_layer'] = self.back_head_layer.as_dict(skip=skip) this_dict['constrain_heads'] = self.constrain_heads this_dict['conformal_roughness'] = self.conformal_roughness diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py index a17ae0b0..31e9d47e 100644 --- a/tests/sample/assemblies/test_bilayer.py +++ b/tests/sample/assemblies/test_bilayer.py @@ -57,7 +57,7 @@ def test_layer_structure(self): def test_custom_layers(self): """Test creation with custom head/tail layers.""" d2o = Material(sld=6.36, isld=0, name='D2O') - air = Material(sld=0, isld=0, name='Air') + air_matched_water = Material(sld=0, isld=0, name='Air Matched Water') front_head = LayerAreaPerMolecule( molecular_formula='C10H18NO8P', @@ -71,7 +71,7 @@ def test_custom_layers(self): tail = LayerAreaPerMolecule( molecular_formula='C32D64', thickness=18.0, - solvent=air, + solvent=air_matched_water, solvent_fraction=0.0, area_per_molecule=50.0, roughness=2.0, @@ -89,7 +89,7 @@ def test_custom_layers(self): p = Bilayer( front_head_layer=front_head, - tail_layer=tail, + front_tail_layer=tail, back_head_layer=back_head, name='Custom Bilayer', ) @@ -275,7 +275,7 @@ def test_dict_round_trip(): """Test serialization/deserialization round trip.""" # When d2o = Material(sld=6.36, isld=0, name='D2O') - air = Material(sld=0, isld=0, name='Air') + air_matched_water = Material(sld=0, isld=0, name='Air Matched Water') front_head = LayerAreaPerMolecule( molecular_formula='C10H18NO8P', @@ -289,7 +289,7 @@ def test_dict_round_trip(): tail = LayerAreaPerMolecule( molecular_formula='C32D64', thickness=18.0, - solvent=air, + solvent=air_matched_water, solvent_fraction=0.0, area_per_molecule=50.0, roughness=2.0, @@ -307,7 +307,7 @@ def test_dict_round_trip(): p = Bilayer( front_head_layer=front_head, - tail_layer=tail, + front_tail_layer=tail, back_head_layer=back_head, constrain_heads=False, conformal_roughness=False, From 4333c43226ba7a2e9b2c351b53450a9d7217f14b Mon Sep 17 00:00:00 2001 From: rozyczko Date: Mon, 2 Feb 2026 15:16:05 +0100 Subject: [PATCH 4/9] Updated notebook after PR --- docs/src/tutorials/simulation/bilayer.ipynb | 332 ++++++++++++++------ tests/sample/assemblies/test_bilayer.py | 8 + 2 files changed, 239 insertions(+), 101 deletions(-) diff --git a/docs/src/tutorials/simulation/bilayer.ipynb b/docs/src/tutorials/simulation/bilayer.ipynb index 5e8d3a0f..77c55ff0 100644 --- a/docs/src/tutorials/simulation/bilayer.ipynb +++ b/docs/src/tutorials/simulation/bilayer.ipynb @@ -145,7 +145,7 @@ "# Create the bilayer with default constraints\n", "bilayer = Bilayer(\n", " front_head_layer=head_layer,\n", - " front_front_front_tail_layer=front_tail_layer,\n", + " front_tail_layer=front_tail_layer,\n", " constrain_heads=True, # Head layers share thickness and area per molecule\n", " conformal_roughness=True, # All layers share the same roughness\n", " name='DPPC Bilayer'\n", @@ -171,16 +171,9 @@ "metadata": {}, "outputs": [], "source": [ - "# The bilayer has 4 layers\n", - "print(f'Number of layers: {len(bilayer.layers)}')\n", - "print()\n", - "\n", - "# Access individual layers\n", - "print('Layer structure (from front to back):')\n", - "print(f' 1. Front Head: {bilayer.front_head_layer.name}')\n", - "print(f' 2. Front Tail: {bilayer.front_tail_layer.name}')\n", - "print(f' 3. Back Tail: {bilayer.back_tail_layer.name}')\n", - "print(f' 4. Back Head: {bilayer.back_head_layer.name}')" + "# The bilayer has 4 layers: front_head, front_tail, back_tail, back_head\n", + "# All structural parameters are automatically constrained\n", + "print(f'Bilayer layers: {len(bilayer.layers)}')" ] }, { @@ -200,19 +193,11 @@ "metadata": {}, "outputs": [], "source": [ - "# Check tail layer constraints - both tails should have the same parameters\n", - "print('Tail layer parameters:')\n", - "print(f' Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", - "print(f' Back tail thickness: {bilayer.back_tail_layer.thickness.value:.2f} Å')\n", - "print(f' Front tail APM: {bilayer.front_tail_layer.area_per_molecule:.2f} Ų')\n", - "print(f' Back tail APM: {bilayer.back_tail_layer.area_per_molecule:.2f} Ų')\n", - "print()\n", - "\n", - "# Change front tail thickness - back should follow\n", - "bilayer.front_tail_layer.thickness.value = 18.0\n", - "print('After changing front tail thickness to 18 Å:')\n", - "print(f' Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", - "print(f' Back tail thickness: {bilayer.back_tail_layer.thickness.value:.2f} Å')" + "# Access key structural parameters\n", + "print(f'Head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n", + "print(f'Tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", + "print(f'Area per molecule: {bilayer.front_head_layer.area_per_molecule:.2f} Ų')\n", + "print(f'Roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')" ] }, { @@ -222,23 +207,10 @@ "metadata": {}, "outputs": [], "source": [ - "# Check head layer constraints\n", - "print('Head layer parameters:')\n", - "print(f' Front head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n", - "print(f' Back head thickness: {bilayer.back_head_layer.thickness.value:.2f} Å')\n", - "print()\n", - "\n", - "# Hydration (solvent fraction) remains independent\n", - "print('Head layer hydration (independent):')\n", - "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", - "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')\n", - "\n", - "# Change back head hydration - front should NOT change\n", - "bilayer.back_head_layer.solvent_fraction = 0.5\n", - "print()\n", - "print('After changing back head solvent fraction to 0.5:')\n", - "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", - "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')" + "# Head layers share thickness and area per molecule\n", + "# But hydration (solvent fraction) remains independent for each side\n", + "print(f'Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", + "print(f'Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')" ] }, { @@ -248,21 +220,9 @@ "metadata": {}, "outputs": [], "source": [ - "# Check conformal roughness\n", - "print('Roughness values (conformal):')\n", - "print(f' Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", - "print(f' Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back tail roughness: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back head roughness: {bilayer.back_head_layer.roughness.value:.2f} Å')\n", - "print()\n", - "\n", - "# Change front head roughness - all should follow\n", - "bilayer.front_head_layer.roughness.value = 4.0\n", - "print('After changing front head roughness to 4.0 Å:')\n", - "print(f' Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", - "print(f' Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back tail roughness: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back head roughness: {bilayer.back_head_layer.roughness.value:.2f} Å')" + "# Conformal roughness: all layers share the same roughness value\n", + "# (controlled by the front head layer)\n", + "print(f'Conformal roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')" ] }, { @@ -287,24 +247,22 @@ "bilayer.front_tail_layer.thickness.value = 16.0\n", "bilayer.back_head_layer.solvent_fraction = 0.3\n", "\n", - "# Create superphase (air) and subphase (silicon substrate with oxide layer)\n", - "vacuum = Material(sld=0, isld=0, name='Vacuum')\n", - "superphase = Layer(material=vacuum, thickness=0, roughness=0, name='Vacuum Superphase')\n", - "\n", + "# Create substrate layers (silicon with oxide layer)\n", "si = Material(sld=2.047, isld=0, name='Si')\n", "sio2 = Material(sld=3.47, isld=0, name='SiO2')\n", - "si_layer = Layer(material=si, thickness=0, roughness=0, name='Si')\n", + "si_layer = Layer(material=si, thickness=0, roughness=0, name='Si Substrate')\n", "sio2_layer = Layer(material=sio2, thickness=15, roughness=3, name='SiO2')\n", "\n", - "# D2O subphase\n", - "d2o_layer = Layer(material=d2o, thickness=0, roughness=3, name='D2O Subphase')\n", + "# D2O subphase (bulk water)\n", + "d2o_subphase = Layer(material=d2o, thickness=0, roughness=3, name='D2O Bulk')\n", "\n", - "# Create sample: superphase -> bilayer -> SiO2 -> Si\n", + "# Create sample structure: Si | SiO2 | head | tail | tail | head | D2O\n", + "# In easyreflectometry convention: superphase (Si) -> layers -> subphase (D2O)\n", "sample = Sample(\n", - " Multilayer(superphase, name='Superphase'),\n", + " Multilayer(si_layer, name='Si Substrate'),\n", + " Multilayer(sio2_layer, name='SiO2'),\n", " bilayer,\n", - " Multilayer(d2o_layer, name='D2O'),\n", - " Multilayer([sio2_layer, si_layer], name='Substrate'),\n", + " Multilayer(d2o_subphase, name='D2O Subphase'),\n", " name='Bilayer on Si/SiO2'\n", ")\n", "\n", @@ -404,20 +362,15 @@ "metadata": {}, "outputs": [], "source": [ - "# Disable conformal roughness\n", + "# Constraints can be modified at runtime\n", + "# Disable conformal roughness to set different values for each layer\n", "bilayer.conformal_roughness = False\n", "\n", - "# Now we can set different roughness values for each layer\n", + "# Set varying roughness across the bilayer\n", "bilayer.front_head_layer.roughness.value = 2.0\n", "bilayer.front_tail_layer.roughness.value = 1.5\n", "bilayer.back_tail_layer.roughness.value = 1.5\n", - "bilayer.back_head_layer.roughness.value = 4.0\n", - "\n", - "print('Individual roughness values after disabling conformal roughness:')\n", - "print(f' Front head: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", - "print(f' Front tail: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back tail: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", - "print(f' Back head: {bilayer.back_head_layer.roughness.value:.2f} Å')" + "bilayer.back_head_layer.roughness.value = 4.0" ] }, { @@ -463,33 +416,207 @@ "metadata": {}, "outputs": [], "source": [ - "# Create a bilayer with asymmetric hydration\n", - "# This is common in supported bilayers where one side is in contact with a substrate\n", + "# Asymmetric hydration: different solvent exposure on each side\n", + "# Set up symmetric hydration first\n", + "bilayer.front_head_layer.solvent_fraction = 0.3\n", + "bilayer.back_head_layer.solvent_fraction = 0.3\n", "\n", - "# Low hydration on front (substrate side)\n", - "bilayer.front_head_layer.solvent_fraction = 0.1\n", + "# Get symmetric SLD profile\n", + "z_sym, sld_sym = model.interface().sld_profile(model.unique_name)\n", "\n", - "# Higher hydration on back (solution side)\n", - "bilayer.back_head_layer.solvent_fraction = 0.4\n", + "# Now set asymmetric hydration (common in supported bilayers)\n", + "bilayer.front_head_layer.solvent_fraction = 0.1 # Substrate side - less hydrated\n", + "bilayer.back_head_layer.solvent_fraction = 0.4 # Solution side - more hydrated\n", "\n", - "print('Asymmetric hydration:')\n", - "print(f' Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", - "print(f' Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')\n", + "# Get asymmetric SLD profile\n", + "z_asym, sld_asym = model.interface().sld_profile(model.unique_name)\n", + "\n", + "# Plot comparison\n", + "plt.figure(figsize=(10, 6))\n", + "plt.plot(z_sym, sld_sym, 'b-', linewidth=2, label='Symmetric (0.3/0.3)')\n", + "plt.plot(z_asym, sld_asym, 'r--', linewidth=2, label='Asymmetric (0.1/0.4)')\n", + "plt.xlabel('Distance from interface (Å)')\n", + "plt.ylabel('SLD (10⁻⁶ Å⁻²)')\n", + "plt.title('Effect of Asymmetric Hydration on SLD Profile')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "df3d351f", + "metadata": {}, + "source": [ + "## Multiple Contrast Analysis\n", "\n", - "# Calculate new reflectometry\n", - "reflectivity_asymmetric = model.interface().reflectity_profile(q, model.unique_name)\n", + "The most common use case for bilayer models is fitting multiple contrast data - measuring the same sample with different isotopic compositions of the solvent (e.g., D₂O, H₂O, or mixtures).\n", "\n", - "# Reset to symmetric for comparison\n", + "The `Bilayer` assembly provides a `constrain_multiple_contrast` method to link structural parameters across different contrast measurements while allowing contrast-specific parameters (like solvent fraction) to vary independently." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2cb3a668", + "metadata": {}, + "outputs": [], + "source": [ + "# Reset bilayer to symmetric hydration\n", "bilayer.front_head_layer.solvent_fraction = 0.3\n", "bilayer.back_head_layer.solvent_fraction = 0.3\n", - "reflectivity_symmetric = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "# Create H2O-matched materials for second contrast\n", + "h2o = Material(sld=-0.56, isld=0, name='H2O')\n", + "air_matched_h2o = Material(sld=0, isld=0, name='Air Matched H2O')\n", + "\n", + "# Create head layer for H2O contrast\n", + "head_layer_h2o = LayerAreaPerMolecule(\n", + " molecular_formula='C10H18NO8P',\n", + " thickness=10.0,\n", + " solvent=h2o,\n", + " solvent_fraction=0.3,\n", + " area_per_molecule=48.2,\n", + " roughness=3.0,\n", + " name='DPPC Head H2O'\n", + ")\n", + "\n", + "# Create tail layer for H2O contrast (hydrogenous tails)\n", + "tail_layer_h2o = LayerAreaPerMolecule(\n", + " molecular_formula='C32H64',\n", + " thickness=16.0,\n", + " solvent=air_matched_h2o,\n", + " solvent_fraction=0.0,\n", + " area_per_molecule=48.2,\n", + " roughness=3.0,\n", + " name='DPPC Tail H2O'\n", + ")\n", + "\n", + "# Create H2O bilayer\n", + "bilayer_h2o = Bilayer(\n", + " front_head_layer=head_layer_h2o,\n", + " front_tail_layer=tail_layer_h2o,\n", + " constrain_heads=True,\n", + " conformal_roughness=True,\n", + " name='DPPC Bilayer H2O'\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b62b234", + "metadata": {}, + "outputs": [], + "source": [ + "# Constrain structural parameters between contrasts\n", + "# Link thicknesses and areas per molecule, but NOT solvent fractions\n", + "bilayer_h2o.constrain_multiple_contrast(\n", + " bilayer,\n", + " front_head_thickness=True,\n", + " back_head_thickness=True,\n", + " tail_thickness=True,\n", + " front_head_area_per_molecule=True,\n", + " back_head_area_per_molecule=True,\n", + " tail_area_per_molecule=True,\n", + " front_head_fraction=False, # Hydration can differ between contrasts\n", + " back_head_fraction=False,\n", + " tail_fraction=False,\n", + ")\n", + "\n", + "print('Structural parameters are now constrained between D2O and H2O contrasts')\n", + "print(f'D2O head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n", + "print(f'H2O head thickness: {bilayer_h2o.front_head_layer.thickness.value:.2f} Å')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a962d9ce", + "metadata": {}, + "outputs": [], + "source": [ + "# Create complete sample for H2O contrast\n", + "h2o_subphase = Layer(material=h2o, thickness=0, roughness=3, name='H2O Subphase')\n", + "\n", + "sample_h2o = Sample(\n", + " Multilayer(si_layer, name='Si Substrate'),\n", + " Multilayer(sio2_layer, name='SiO2'),\n", + " bilayer_h2o,\n", + " Multilayer(h2o_subphase, name='H2O Subphase'),\n", + " name='Bilayer on Si/SiO2 in H2O'\n", + ")\n", + "\n", + "# Create model for H2O contrast\n", + "model_h2o = Model(\n", + " sample=sample_h2o,\n", + " scale=1.0,\n", + " background=1e-7,\n", + " resolution_function=PercentageFwhm(5),\n", + " name='Bilayer Model H2O'\n", + ")\n", + "model_h2o.interface = interface" + ] + }, + { + "cell_type": "markdown", + "id": "75f50b75", + "metadata": {}, + "source": [ + "### SLD Profiles: D₂O vs H₂O Contrast\n", + "\n", + "The SLD profiles clearly show the difference in contrast between D₂O and H₂O measurements." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1676ddab", + "metadata": {}, + "outputs": [], + "source": [ + "# Get SLD profiles for both contrasts\n", + "z_d2o, sld_d2o = model.interface().sld_profile(model.unique_name)\n", + "z_h2o, sld_h2o = model_h2o.interface().sld_profile(model_h2o.unique_name)\n", "\n", "plt.figure(figsize=(10, 6))\n", - "plt.semilogy(q, reflectivity_symmetric, 'b-', linewidth=2, label='Symmetric hydration (0.3/0.3)')\n", - "plt.semilogy(q, reflectivity_asymmetric, 'r--', linewidth=2, label='Asymmetric hydration (0.1/0.4)')\n", + "plt.plot(z_d2o, sld_d2o, 'b-', linewidth=2, label='D₂O contrast')\n", + "plt.plot(z_h2o, sld_h2o, 'r-', linewidth=2, label='H₂O contrast')\n", + "plt.xlabel('Distance from interface (Å)')\n", + "plt.ylabel('SLD (10⁻⁶ Å⁻²)')\n", + "plt.title('SLD Profiles: Multiple Contrast Comparison')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "5da06bc3", + "metadata": {}, + "source": [ + "### Reflectivity Curves: D₂O vs H₂O Contrast\n", + "\n", + "The reflectivity curves show how the different contrasts provide complementary structural information." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b5688bc", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate reflectivity for both contrasts\n", + "reflectivity_d2o = model.interface().reflectity_profile(q, model.unique_name)\n", + "reflectivity_h2o = model_h2o.interface().reflectity_profile(q, model_h2o.unique_name)\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.semilogy(q, reflectivity_d2o, 'b-', linewidth=2, label='D₂O contrast')\n", + "plt.semilogy(q, reflectivity_h2o, 'r-', linewidth=2, label='H₂O contrast')\n", "plt.xlabel('Q (Å⁻¹)')\n", "plt.ylabel('Reflectivity')\n", - "plt.title('Effect of Asymmetric Hydration on Reflectometry')\n", + "plt.title('Reflectivity: Multiple Contrast Comparison')\n", "plt.legend()\n", "plt.grid(True, alpha=0.3)\n", "plt.show()" @@ -497,24 +624,27 @@ }, { "cell_type": "markdown", - "id": "df3d351f", + "id": "624b1d60", "metadata": {}, "source": [ "## Summary\n", "\n", - "In this tutorial, we have explored the `Bilayer` assembly in `easyreflectometry`:\n", + "In this tutorial, we explored the `Bilayer` assembly in `easyreflectometry`:\n", "\n", "1. **Creating a bilayer**: Using `LayerAreaPerMolecule` components for head and tail layers\n", - "2. **Understanding constraints**: \n", + "2. **Built-in constraints**: \n", " - Tail layers share all structural parameters\n", " - Head layers share thickness and area per molecule (hydration is independent)\n", " - Conformal roughness applies to all layers by default\n", - "3. **Building a sample**: Combining the bilayer with sub- and super-phases\n", - "4. **Simulating reflectometry**: Using the calculator interface to generate reflectivity profiles\n", - "5. **Modifying constraints**: Disabling conformal roughness for more complex models\n", - "6. **Asymmetric hydration**: Modeling supported bilayers with different solvent exposure on each side\n", - "\n", - "The `Bilayer` assembly provides a convenient way to model phospholipid bilayers with physically meaningful constraints, reducing the number of free parameters while maintaining flexibility for complex systems." + "3. **Building a sample**: Combining the bilayer with sub- and super-phases (Si/SiO₂ substrate, water subphase)\n", + "4. **Simulating reflectometry**: Using the calculator interface to generate reflectivity and SLD profiles\n", + "5. **Asymmetric hydration**: Modeling supported bilayers with different solvent exposure on each side\n", + "6. **Multiple contrast analysis**: \n", + " - Creating bilayers with different solvent contrasts (D₂O/H₂O)\n", + " - Using `constrain_multiple_contrast` to link structural parameters across contrasts\n", + " - Visualizing complementary information from different contrasts\n", + "\n", + "The `Bilayer` assembly provides a convenient way to model phospholipid bilayers with physically meaningful constraints, making it ideal for simultaneous fitting of multiple contrast reflectometry data." ] } ], diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py index 31e9d47e..e565f1e4 100644 --- a/tests/sample/assemblies/test_bilayer.py +++ b/tests/sample/assemblies/test_bilayer.py @@ -15,6 +15,14 @@ class TestBilayer: + def setup_method(self): + from easyscience import global_object + # Clear the global object map to prevent name collisions + # Accessing private _clear method as Map doesn't expose a public clear + if hasattr(global_object.map, 'clear'): + global_object.map.clear() + elif hasattr(global_object.map, '_clear'): + global_object.map._clear() def test_default(self): """Test default bilayer creation with expected structure.""" p = Bilayer() From b703e8339d921e01f7f19efb6f384fb3f560beaf Mon Sep 17 00:00:00 2001 From: rozyczko Date: Thu, 12 Feb 2026 14:59:03 +0100 Subject: [PATCH 5/9] functionality review issues addressed --- docs/src/tutorials/simulation/bilayer.ipynb | 96 +++++++++++++-------- 1 file changed, 60 insertions(+), 36 deletions(-) diff --git a/docs/src/tutorials/simulation/bilayer.ipynb b/docs/src/tutorials/simulation/bilayer.ipynb index 77c55ff0..3d7faccd 100644 --- a/docs/src/tutorials/simulation/bilayer.ipynb +++ b/docs/src/tutorials/simulation/bilayer.ipynb @@ -65,7 +65,7 @@ "\n", "We'll create a DPPC (dipalmitoylphosphatidylcholine) bilayer, a common model phospholipid.\n", "\n", - "First, let's define the materials for our solvents - D₂O (heavy water) for the head groups and air for the tail groups." + "First, let's define the solvent material — D₂O (heavy water) — which is used for all layers in the bilayer." ] }, { @@ -75,9 +75,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Define solvent materials\n", - "d2o = Material(sld=6.36, isld=0, name='D2O')\n", - "air_matched_water = Material(sld=0, isld=0, name='Air Matched Water')" + "# Define the solvent material\n", + "d2o = Material(sld=6.36, isld=0, name='D2O')" ] }, { @@ -114,7 +113,7 @@ "front_tail_layer = LayerAreaPerMolecule(\n", " molecular_formula='C32D64',\n", " thickness=16.0,\n", - " solvent=air_matched_water,\n", + " solvent=d2o,\n", " solvent_fraction=0.0, # No solvent in the tail region\n", " area_per_molecule=48.2,\n", " roughness=3.0,\n", @@ -183,7 +182,7 @@ "source": [ "### Verifying Constraints\n", "\n", - "Let's verify that the constraints are working correctly." + "The `Bilayer` assembly provides access to all four sub-layers through the properties `front_head_layer`, `front_tail_layer`, `back_tail_layer`, and `back_head_layer`. The back tail and back head layers are automatically created by the assembly, with their structural parameters constrained to the corresponding front layers. Let's verify that the constraints are working correctly." ] }, { @@ -196,8 +195,17 @@ "# Access key structural parameters\n", "print(f'Head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n", "print(f'Tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n", - "print(f'Area per molecule: {bilayer.front_head_layer.area_per_molecule:.2f} Ų')\n", - "print(f'Roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')" + "print(f'Area per molecule: {bilayer.front_head_layer.area_per_molecule:.2f} Ų')" + ] + }, + { + "cell_type": "markdown", + "id": "219395b2", + "metadata": {}, + "source": [ + "### Independent Head Layer Hydration\n", + "\n", + "While `constrain_heads=True` links the head layers' thicknesses and areas per molecule, the solvent fraction (hydration) of each head layer remains independent. This is physically meaningful — in a supported bilayer, the head group facing the substrate may have different hydration than the one facing the bulk solvent. We can access the front and back head layers through `bilayer.front_head_layer` and `bilayer.back_head_layer`." ] }, { @@ -207,12 +215,28 @@ "metadata": {}, "outputs": [], "source": [ - "# Head layers share thickness and area per molecule\n", - "# But hydration (solvent fraction) remains independent for each side\n", + "# Head layers share thickness and area per molecule via constrain_heads=True,\n", + "# but solvent fraction is independent and can be set separately for each side.\n", + "print(f'Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", + "print(f'Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')\n", + "\n", + "# We can set them independently\n", + "bilayer.back_head_layer.solvent_fraction = 0.5\n", + "print(f'\\nAfter setting back head solvent fraction to 0.5:')\n", "print(f'Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n", "print(f'Back head solvent fraction: {bilayer.back_head_layer.solvent_fraction:.2f}')" ] }, + { + "cell_type": "markdown", + "id": "b278da84", + "metadata": {}, + "source": [ + "### Conformal Roughness\n", + "\n", + "When `conformal_roughness=True`, all four layers in the bilayer share the same roughness value, controlled by the front head layer's roughness parameter. Let's verify this by printing the roughness of each layer." + ] + }, { "cell_type": "code", "execution_count": null, @@ -222,7 +246,10 @@ "source": [ "# Conformal roughness: all layers share the same roughness value\n", "# (controlled by the front head layer)\n", - "print(f'Conformal roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')" + "print(f'Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n", + "print(f'Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n", + "print(f'Back tail roughness: {bilayer.back_tail_layer.roughness.value:.2f} Å')\n", + "print(f'Back head roughness: {bilayer.back_head_layer.roughness.value:.2f} Å')" ] }, { @@ -362,32 +389,25 @@ "metadata": {}, "outputs": [], "source": [ - "# Constraints can be modified at runtime\n", - "# Disable conformal roughness to set different values for each layer\n", + "# First, compute reflectivity with current conformal roughness (3.0 Å)\n", + "reflectivity_conformal = model.interface().reflectity_profile(q, model.unique_name)\n", + "\n", + "# Disable conformal roughness to allow independent roughness per layer\n", "bilayer.conformal_roughness = False\n", "\n", + "# Re-establish calculator bindings after constraint change\n", + "model.generate_bindings()\n", + "\n", "# Set varying roughness across the bilayer\n", "bilayer.front_head_layer.roughness.value = 2.0\n", "bilayer.front_tail_layer.roughness.value = 1.5\n", "bilayer.back_tail_layer.roughness.value = 1.5\n", - "bilayer.back_head_layer.roughness.value = 4.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5fd26d95", - "metadata": {}, - "outputs": [], - "source": [ - "# Compare reflectometry with different roughness configurations\n", - "reflectivity_variable_roughness = model.interface().reflectity_profile(q, model.unique_name)\n", + "bilayer.back_head_layer.roughness.value = 4.0\n", "\n", - "# Reset to conformal roughness\n", - "bilayer.conformal_roughness = True\n", - "bilayer.front_head_layer.roughness.value = 3.0\n", - "reflectivity_conformal = model.interface().reflectity_profile(q, model.unique_name)\n", + "# Compute reflectivity with variable roughness\n", + "reflectivity_variable_roughness = model.interface().reflectity_profile(q, model.unique_name)\n", "\n", + "# Plot comparison\n", "plt.figure(figsize=(10, 6))\n", "plt.semilogy(q, reflectivity_conformal, 'b-', linewidth=2, label='Conformal roughness (3.0 Å)')\n", "plt.semilogy(q, reflectivity_variable_roughness, 'r--', linewidth=2, label='Variable roughness')\n", @@ -396,7 +416,12 @@ "plt.title('Effect of Roughness Configuration on Reflectometry')\n", "plt.legend()\n", "plt.grid(True, alpha=0.3)\n", - "plt.show()" + "plt.show()\n", + "\n", + "# Reset to conformal roughness for subsequent cells\n", + "bilayer.conformal_roughness = True\n", + "bilayer.front_head_layer.roughness.value = 3.0\n", + "model.generate_bindings()" ] }, { @@ -466,11 +491,10 @@ "bilayer.front_head_layer.solvent_fraction = 0.3\n", "bilayer.back_head_layer.solvent_fraction = 0.3\n", "\n", - "# Create H2O-matched materials for second contrast\n", + "# Create H2O material for second contrast\n", "h2o = Material(sld=-0.56, isld=0, name='H2O')\n", - "air_matched_h2o = Material(sld=0, isld=0, name='Air Matched H2O')\n", "\n", - "# Create head layer for H2O contrast\n", + "# Create head layer for H2O contrast (same lipid, different solvent)\n", "head_layer_h2o = LayerAreaPerMolecule(\n", " molecular_formula='C10H18NO8P',\n", " thickness=10.0,\n", @@ -481,11 +505,11 @@ " name='DPPC Head H2O'\n", ")\n", "\n", - "# Create tail layer for H2O contrast (hydrogenous tails)\n", + "# Create tail layer for H2O contrast (same deuterated lipid, different solvent)\n", "tail_layer_h2o = LayerAreaPerMolecule(\n", - " molecular_formula='C32H64',\n", + " molecular_formula='C32D64',\n", " thickness=16.0,\n", - " solvent=air_matched_h2o,\n", + " solvent=h2o,\n", " solvent_fraction=0.0,\n", " area_per_molecule=48.2,\n", " roughness=3.0,\n", From 7217924ec1403542f8a938555fb978c24e6f1e12 Mon Sep 17 00:00:00 2001 From: rozyczko Date: Thu, 12 Feb 2026 15:27:01 +0100 Subject: [PATCH 6/9] more PR review comments addressed --- .../sample/assemblies/bilayer.py | 297 +++++++++++------- tests/sample/assemblies/test_bilayer.py | 4 +- 2 files changed, 186 insertions(+), 115 deletions(-) diff --git a/src/easyreflectometry/sample/assemblies/bilayer.py b/src/easyreflectometry/sample/assemblies/bilayer.py index 797ca792..21c428f6 100644 --- a/src/easyreflectometry/sample/assemblies/bilayer.py +++ b/src/easyreflectometry/sample/assemblies/bilayer.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import Optional +from typing import Any from easyscience import global_object from easyscience.variable import Parameter @@ -10,6 +10,28 @@ from ..elements.materials.material import Material from .base_assembly import BaseAssembly +DEFAULTS = { + 'head': { + 'molecular_formula': 'C10H18NO8P', + 'thickness': 10.0, + 'solvent_fraction': 0.2, + 'area_per_molecule': 48.2, + 'roughness': 3.0, + }, + 'tail': { + 'molecular_formula': 'C32D64', + 'thickness': 16.0, + 'solvent_fraction': 0.0, + 'area_per_molecule': 48.2, + 'roughness': 3.0, + }, + 'solvent': { + 'sld': 6.36, + 'isld': 0, + 'name': 'D2O', + }, +} + class Bilayer(BaseAssembly): """A lipid bilayer consisting of two surfactant layers where one is inverted. @@ -17,9 +39,11 @@ class Bilayer(BaseAssembly): The bilayer structure is: Front Head - Front Tail - Back Tail - Back Head This assembly comes pre-populated with physically meaningful constraints: - - Both tail layers share the same structural parameters (thickness, area per molecule) - - Head layers share thickness and area per molecule (different hydration/solvent_fraction allowed) - - A single roughness parameter applies to all interfaces (conformal roughness) + - Both tail layers are constrained to share the same structural parameters + (thickness, area per molecule, and solvent fraction). + - Head layers are constrained to share thickness and area per molecule, + while solvent fraction (hydration) remains independent on each side. + - A single roughness parameter applies to all interfaces (conformal roughness). More information about the usage of this assembly is available in the `bilayer documentation`_ @@ -29,118 +53,64 @@ class Bilayer(BaseAssembly): def __init__( self, - front_head_layer: Optional[LayerAreaPerMolecule] = None, - front_tail_layer: Optional[LayerAreaPerMolecule] = None, - back_head_layer: Optional[LayerAreaPerMolecule] = None, + front_head_layer: LayerAreaPerMolecule | None = None, + front_tail_layer: LayerAreaPerMolecule | None = None, + back_head_layer: LayerAreaPerMolecule | None = None, name: str = 'EasyBilayer', - unique_name: Optional[str] = None, + unique_name: str | None = None, constrain_heads: bool = True, conformal_roughness: bool = True, - interface=None, + interface: Any = None, ): """Constructor. :param front_head_layer: Layer representing the front head part of the bilayer. - :param front_tail_layer: Layer representing the tail part of the bilayer. A second tail - layer is created internally with parameters constrained to this one. + :param front_tail_layer: Layer representing the front tail part of the bilayer. + A back tail layer is created internally with its thickness, area per molecule, + and solvent fraction constrained to match this layer. :param back_head_layer: Layer representing the back head part of the bilayer. :param name: Name for bilayer, defaults to 'EasyBilayer'. - :param constrain_heads: Constrain head layer thickness and area per molecule, defaults to `True`. - :param conformal_roughness: Constrain roughness to be the same for all layers, defaults to `True`. + :param unique_name: Unique name for internal object tracking, defaults to `None`. + :param constrain_heads: When `True`, the back head layer thickness and area per + molecule are constrained to match the front head layer. Solvent fraction + (hydration) remains independent on each side. Defaults to `True`. + :param conformal_roughness: When `True`, all four layer interfaces share + the same roughness value, controlled by the front head layer. Defaults to `True`. :param interface: Calculator interface, defaults to `None`. """ # Generate unique name for nested objects if unique_name is None: unique_name = global_object.generate_unique_name(self.__class__.__name__) - # Create default front head layer if not provided + # Create default layers if not provided if front_head_layer is None: - d2o_front = Material( - sld=6.36, - isld=0, - name='D2O', - unique_name=unique_name + '_MaterialFrontHead', - interface=interface, - ) - front_head_layer = LayerAreaPerMolecule( - molecular_formula='C10H18NO8P', - thickness=10.0, - solvent=d2o_front, - solvent_fraction=0.2, - area_per_molecule=48.2, - roughness=3.0, - name='DPPC Head Front', - unique_name=unique_name + '_LayerAreaPerMoleculeFrontHead', + front_head_layer = self._create_default_head_layer( + unique_name=unique_name, + name_suffix='Front', interface=interface, ) - # Create default front tail layer if not provided if front_tail_layer is None: - d2o_tail = Material( - sld=6.36, - isld=0, - name='D2O', - unique_name=unique_name + '_MaterialTail', - interface=interface, - ) - front_tail_layer = LayerAreaPerMolecule( - molecular_formula='C32D64', - thickness=16.0, - solvent=d2o_tail, - solvent_fraction=0.0, - area_per_molecule=48.2, - roughness=3.0, - name='DPPC Tail', - unique_name=unique_name + '_LayerAreaPerMoleculeTail', + front_tail_layer = self._create_default_tail_layer( + unique_name=unique_name, interface=interface, ) - # Create second tail layer with same parameters as the first - # This will be constrained to the first tail layer - d2o_back_tail = Material( - sld=6.36, - isld=0, - name='D2O', - unique_name=unique_name + '_MaterialBackTail', - interface=interface, - ) - back_tail_layer = LayerAreaPerMolecule( - molecular_formula=front_tail_layer.molecular_formula, - thickness=front_tail_layer.thickness.value, - solvent=d2o_back_tail, - solvent_fraction=front_tail_layer.solvent_fraction, - area_per_molecule=front_tail_layer.area_per_molecule, - roughness=front_tail_layer.roughness.value, - name=front_tail_layer.name + ' Back', - unique_name=unique_name + '_LayerAreaPerMoleculeBackTail', + # Create back tail layer with initial values copied from the front tail. + # Its parameters will be constrained to the front tail after construction. + back_tail_layer = self._create_back_tail_layer( + front_tail_layer=front_tail_layer, + unique_name=unique_name, interface=interface, ) - # Create default back head layer if not provided if back_head_layer is None: - d2o_back = Material( - sld=6.36, - isld=0, - name='D2O', - unique_name=unique_name + '_MaterialBackHead', - interface=interface, - ) - back_head_layer = LayerAreaPerMolecule( - molecular_formula='C10H18NO8P', - thickness=10.0, - solvent=d2o_back, - solvent_fraction=0.2, - area_per_molecule=48.2, - roughness=3.0, - name='DPPC Head Back', - unique_name=unique_name + '_LayerAreaPerMoleculeBackHead', + back_head_layer = self._create_default_head_layer( + unique_name=unique_name, + name_suffix='Back', interface=interface, ) - # Store the front tail layer reference for constraint setup - self._front_tail_layer = front_tail_layer - self._back_tail_layer = back_tail_layer - # Create layer collection: front_head, front_tail, back_tail, back_head bilayer_layers = LayerCollection( front_head_layer, @@ -176,10 +146,108 @@ def __init__( if conformal_roughness: self.conformal_roughness = True + @staticmethod + def _create_default_head_layer( + unique_name: str, + name_suffix: str, + interface: Any = None, + ) -> LayerAreaPerMolecule: + """Create a default head layer with DPPC head group parameters. + + :param unique_name: Base unique name for internal object tracking. + :param name_suffix: Suffix for layer name ('Front' or 'Back'). + :param interface: Calculator interface, defaults to `None`. + :return: A new LayerAreaPerMolecule for the head group. + """ + solvent = Material( + sld=DEFAULTS['solvent']['sld'], + isld=DEFAULTS['solvent']['isld'], + name=DEFAULTS['solvent']['name'], + unique_name=unique_name + f'_Material{name_suffix}Head', + interface=interface, + ) + return LayerAreaPerMolecule( + molecular_formula=DEFAULTS['head']['molecular_formula'], + thickness=DEFAULTS['head']['thickness'], + solvent=solvent, + solvent_fraction=DEFAULTS['head']['solvent_fraction'], + area_per_molecule=DEFAULTS['head']['area_per_molecule'], + roughness=DEFAULTS['head']['roughness'], + name=f'DPPC Head {name_suffix}', + unique_name=unique_name + f'_LayerAreaPerMolecule{name_suffix}Head', + interface=interface, + ) + + @staticmethod + def _create_default_tail_layer( + unique_name: str, + interface: Any = None, + ) -> LayerAreaPerMolecule: + """Create a default tail layer with DPPC tail group parameters. + + :param unique_name: Base unique name for internal object tracking. + :param interface: Calculator interface, defaults to `None`. + :return: A new LayerAreaPerMolecule for the tail group. + """ + solvent = Material( + sld=DEFAULTS['solvent']['sld'], + isld=DEFAULTS['solvent']['isld'], + name=DEFAULTS['solvent']['name'], + unique_name=unique_name + '_MaterialTail', + interface=interface, + ) + return LayerAreaPerMolecule( + molecular_formula=DEFAULTS['tail']['molecular_formula'], + thickness=DEFAULTS['tail']['thickness'], + solvent=solvent, + solvent_fraction=DEFAULTS['tail']['solvent_fraction'], + area_per_molecule=DEFAULTS['tail']['area_per_molecule'], + roughness=DEFAULTS['tail']['roughness'], + name='DPPC Tail', + unique_name=unique_name + '_LayerAreaPerMoleculeTail', + interface=interface, + ) + + @staticmethod + def _create_back_tail_layer( + front_tail_layer: LayerAreaPerMolecule, + unique_name: str, + interface: Any = None, + ) -> LayerAreaPerMolecule: + """Create a back tail layer with initial values copied from the front tail layer. + + :param front_tail_layer: The front tail layer to copy initial values from. + :param unique_name: Base unique name for internal object tracking. + :param interface: Calculator interface, defaults to `None`. + :return: A new LayerAreaPerMolecule for the back tail. + """ + solvent = Material( + sld=DEFAULTS['solvent']['sld'], + isld=DEFAULTS['solvent']['isld'], + name=DEFAULTS['solvent']['name'], + unique_name=unique_name + '_MaterialBackTail', + interface=interface, + ) + return LayerAreaPerMolecule( + molecular_formula=front_tail_layer.molecular_formula, + thickness=front_tail_layer.thickness.value, + solvent=solvent, + solvent_fraction=front_tail_layer.solvent_fraction, + area_per_molecule=front_tail_layer.area_per_molecule, + roughness=front_tail_layer.roughness.value, + name=front_tail_layer.name + ' Back', + unique_name=unique_name + '_LayerAreaPerMoleculeBackTail', + interface=interface, + ) + def _setup_tail_constraints(self) -> None: - """Setup constraints so back tail layer parameters depend on front tail layer.""" - front_tail = self._front_tail_layer - back_tail = self._back_tail_layer + """Setup constraints so back tail layer parameters depend on front tail layer. + + Constrains thickness, area per molecule, and solvent fraction of the + back tail layer to match the front tail layer. + """ + front_tail = self.front_tail_layer + back_tail = self.back_tail_layer # Constrain thickness back_tail.thickness.make_dependent_on( @@ -188,21 +256,21 @@ def _setup_tail_constraints(self) -> None: ) # Constrain area per molecule - back_tail._area_per_molecule.make_dependent_on( + back_tail.area_per_molecule_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': front_tail._area_per_molecule}, + dependency_map={'a': front_tail.area_per_molecule_parameter}, ) # Constrain solvent fraction - back_tail.material._fraction.make_dependent_on( + back_tail.solvent_fraction_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': front_tail.material._fraction}, + dependency_map={'a': front_tail.solvent_fraction_parameter}, ) self._tail_constraints_setup = True @property - def front_head_layer(self) -> Optional[LayerAreaPerMolecule]: + def front_head_layer(self) -> LayerAreaPerMolecule: """Get the front head layer of the bilayer.""" return self.layers[0] @@ -212,17 +280,17 @@ def front_head_layer(self, layer: LayerAreaPerMolecule) -> None: self.layers[0] = layer @property - def front_tail_layer(self) -> Optional[LayerAreaPerMolecule]: + def front_tail_layer(self) -> LayerAreaPerMolecule: """Get the front tail layer of the bilayer.""" return self.layers[1] @property - def back_tail_layer(self) -> Optional[LayerAreaPerMolecule]: + def back_tail_layer(self) -> LayerAreaPerMolecule: """Get the back tail layer of the bilayer.""" return self.layers[2] @property - def back_head_layer(self) -> Optional[LayerAreaPerMolecule]: + def back_head_layer(self) -> LayerAreaPerMolecule: """Get the back head layer of the bilayer.""" return self.layers[3] @@ -264,15 +332,15 @@ def _enable_head_constraints(self) -> None: ) # Constrain area per molecule - back_head._area_per_molecule.make_dependent_on( + back_head.area_per_molecule_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': front_head._area_per_molecule}, + dependency_map={'a': front_head.area_per_molecule_parameter}, ) def _disable_head_constraints(self) -> None: """Disable head layer constraints.""" self.back_head_layer.thickness.make_independent() - self.back_head_layer._area_per_molecule.make_independent() + self.back_head_layer.area_per_molecule_parameter.make_independent() @property def conformal_roughness(self) -> bool: @@ -324,6 +392,9 @@ def constrain_multiple_contrast( ) -> None: """Constrain structural parameters between bilayer objects. + Makes this bilayer's parameters dependent on another_contrast's parameters, + so that changes to another_contrast propagate to this bilayer. + :param another_contrast: The bilayer to constrain to. :param front_head_thickness: Constrain front head thickness. :param back_head_thickness: Constrain back head thickness. @@ -354,39 +425,39 @@ def constrain_multiple_contrast( ) if front_head_area_per_molecule: - self.front_head_layer._area_per_molecule.make_dependent_on( + self.front_head_layer.area_per_molecule_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.front_head_layer._area_per_molecule}, + dependency_map={'a': another_contrast.front_head_layer.area_per_molecule_parameter}, ) if back_head_area_per_molecule: - self.back_head_layer._area_per_molecule.make_dependent_on( + self.back_head_layer.area_per_molecule_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.back_head_layer._area_per_molecule}, + dependency_map={'a': another_contrast.back_head_layer.area_per_molecule_parameter}, ) if tail_area_per_molecule: - self.front_tail_layer._area_per_molecule.make_dependent_on( + self.front_tail_layer.area_per_molecule_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.front_tail_layer._area_per_molecule}, + dependency_map={'a': another_contrast.front_tail_layer.area_per_molecule_parameter}, ) if front_head_fraction: - self.front_head_layer.material._fraction.make_dependent_on( + self.front_head_layer.solvent_fraction_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.front_head_layer.material._fraction}, + dependency_map={'a': another_contrast.front_head_layer.solvent_fraction_parameter}, ) if back_head_fraction: - self.back_head_layer.material._fraction.make_dependent_on( + self.back_head_layer.solvent_fraction_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.back_head_layer.material._fraction}, + dependency_map={'a': another_contrast.back_head_layer.solvent_fraction_parameter}, ) if tail_fraction: - self.front_tail_layer.material._fraction.make_dependent_on( + self.front_tail_layer.solvent_fraction_parameter.make_dependent_on( dependency_expression='a', - dependency_map={'a': another_contrast.front_tail_layer.material._fraction}, + dependency_map={'a': another_contrast.front_tail_layer.solvent_fraction_parameter}, ) @property @@ -403,7 +474,7 @@ def _dict_repr(self) -> dict: } } - def as_dict(self, skip: Optional[list[str]] = None) -> dict: + def as_dict(self, skip: list[str] | None = None) -> dict: """Produce a cleaned dict using a custom as_dict method. The resulting dict matches the parameters in __init__ diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py index e565f1e4..1c577c94 100644 --- a/tests/sample/assemblies/test_bilayer.py +++ b/tests/sample/assemblies/test_bilayer.py @@ -123,7 +123,7 @@ def test_tail_layers_linked(self): assert p.back_tail_layer.thickness.value == 20.0 # Change front tail area per molecule - back tail should follow - p.front_tail_layer._area_per_molecule.value = 55.0 + p.front_tail_layer.area_per_molecule = 55.0 assert p.front_tail_layer.area_per_molecule == 55.0 assert p.back_tail_layer.area_per_molecule == 55.0 @@ -137,7 +137,7 @@ def test_constrain_heads_enabled(self): assert p.back_head_layer.thickness.value == 15.0 # Change front head area per molecule - back head should follow - p.front_head_layer._area_per_molecule.value = 60.0 + p.front_head_layer.area_per_molecule = 60.0 assert p.front_head_layer.area_per_molecule == 60.0 assert p.back_head_layer.area_per_molecule == 60.0 From 70066df98ed7f274b3a4da74d20a85b2bc983687 Mon Sep 17 00:00:00 2001 From: Piotr Rozyczko Date: Fri, 13 Feb 2026 19:15:13 +0100 Subject: [PATCH 7/9] Append sample (#289) * fix for appending samples from loaded models * added unit tests * use most recent version of refl1d * update docs to include add_sample_from_orso * temporarily revert to fwhm * proper conversion from stack to layers. * fixed SLD read. * added handling for orso names (model and experiment) * model reindexing issue fixed. Project rst added * path for SLD-less orso files --- docs/src/api/api.rst | 9 + docs/src/api/project.rst | 4 + pyproject.toml | 6 +- .../calculators/calculator_base.py | 2 +- src/easyreflectometry/data/data_store.py | 4 +- src/easyreflectometry/data/measurement.py | 15 +- src/easyreflectometry/fitting.py | 6 +- src/easyreflectometry/model/model.py | 15 + src/easyreflectometry/orso_utils.py | 149 ++++-- src/easyreflectometry/project.py | 231 ++++++++-- .../refl1d/test_refl1d_calculator.py | 4 +- .../calculators/refl1d/test_refl1d_wrapper.py | 6 +- tests/data/test_data_store.py | 19 +- tests/summary/test_summary.py | 4 +- tests/test_data.py | 6 +- tests/test_fitting.py | 36 +- tests/test_measurement_comprehensive.py | 129 +++--- tests/test_orso_utils.py | 131 +++++- tests/test_ort_file.py | 94 ++-- tests/test_project.py | 424 +++++++++++++++++- 20 files changed, 1033 insertions(+), 261 deletions(-) create mode 100644 docs/src/api/project.rst diff --git a/docs/src/api/api.rst b/docs/src/api/api.rst index 6c1d5b3b..ea8e2661 100644 --- a/docs/src/api/api.rst +++ b/docs/src/api/api.rst @@ -19,6 +19,15 @@ Sample is build from assemblies. sample +Project +======= +Project provides a higher-level interface for managing models, experiments, and ORSO import. + +.. toctree:: + :maxdepth: 1 + + project + Assemblies ========== Assemblies are collections of layers that are used to represent a specific physical setup. diff --git a/docs/src/api/project.rst b/docs/src/api/project.rst new file mode 100644 index 00000000..2f8f2932 --- /dev/null +++ b/docs/src/api/project.rst @@ -0,0 +1,4 @@ +.. automodule:: easyreflectometry.project + :members: + :undoc-members: + :show-inheritance: diff --git a/pyproject.toml b/pyproject.toml index 984758ed..630422ca 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,11 +29,11 @@ classifiers = [ requires-python = ">=3.11,<3.13" dependencies = [ - #"easyscience @ git+https://github.com/easyscience/corelib.git@dict_size_changed_bug", - "easyscience", + "easyscience @ git+https://github.com/easyscience/corelib.git@develop", + #"easyscience", "scipp", "refnx", - "refl1d>=1.0.0rc0", + "refl1d>=1.0.0", "orsopy", "svglib<1.6 ; platform_system=='Linux'", "xhtml2pdf", diff --git a/src/easyreflectometry/calculators/calculator_base.py b/src/easyreflectometry/calculators/calculator_base.py index 1be7b568..7d1314cd 100644 --- a/src/easyreflectometry/calculators/calculator_base.py +++ b/src/easyreflectometry/calculators/calculator_base.py @@ -7,7 +7,7 @@ from easyscience.fitting.calculators.interface_factory import ItemContainer from easyscience.io import SerializerComponent -#if TYPE_CHECKING: +# if TYPE_CHECKING: from easyreflectometry.model import Model from easyreflectometry.sample import BaseAssembly from easyreflectometry.sample import Layer diff --git a/src/easyreflectometry/data/data_store.py b/src/easyreflectometry/data/data_store.py index 8dde0bf7..42fd28b0 100644 --- a/src/easyreflectometry/data/data_store.py +++ b/src/easyreflectometry/data/data_store.py @@ -76,7 +76,7 @@ def __init__( y: Optional[Union[np.ndarray, list]] = None, ye: Optional[Union[np.ndarray, list]] = None, xe: Optional[Union[np.ndarray, list]] = None, - model: Optional['Model'] = None, # delay type checking until runtime (quotes) + model: Optional['Model'] = None, # delay type checking until runtime (quotes) x_label: str = 'x', y_label: str = 'y', ): @@ -117,7 +117,7 @@ def __init__( self._color = None @property - def model(self) -> 'Model': # delay type checking until runtime (quotes) + def model(self) -> 'Model': # delay type checking until runtime (quotes) return self._model @model.setter diff --git a/src/easyreflectometry/data/measurement.py b/src/easyreflectometry/data/measurement.py index b24f32a8..df4064b6 100644 --- a/src/easyreflectometry/data/measurement.py +++ b/src/easyreflectometry/data/measurement.py @@ -30,12 +30,25 @@ def load_as_dataset(fname: Union[TextIO, str]) -> DataSet1D: coords_name = 'Qz_' + basename coords_name = list(data_group['coords'].keys())[0] if coords_name not in data_group['coords'] else coords_name data_name = list(data_group['data'].keys())[0] if data_name not in data_group['data'] else data_name - return DataSet1D( + dataset = DataSet1D( x=data_group['coords'][coords_name].values, y=data_group['data'][data_name].values, ye=data_group['data'][data_name].variances, xe=data_group['coords'][coords_name].variances, ) + return dataset + + +def extract_orso_title(data_group: sc.DataGroup, data_name: str) -> str | None: + try: + header = data_group['attrs'][data_name]['orso_header'] + title = header.values.get('data_source', {}).get('experiment', {}).get('title') + except (AttributeError, KeyError, TypeError): + return None + if title is None: + return None + title_str = str(title).strip() + return title_str or None def _load_txt(fname: Union[TextIO, str]) -> sc.DataGroup: diff --git a/src/easyreflectometry/fitting.py b/src/easyreflectometry/fitting.py index 99b5e614..a33e6570 100644 --- a/src/easyreflectometry/fitting.py +++ b/src/easyreflectometry/fitting.py @@ -55,13 +55,13 @@ def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup: variances = data['data'][f'R_{i}'].variances # Find points with non-zero variance - zero_variance_mask = (variances == 0.0) + zero_variance_mask = variances == 0.0 num_zero_variance = np.sum(zero_variance_mask) if num_zero_variance > 0: warnings.warn( - f"Masked {num_zero_variance} data point(s) in reflectivity {i} due to zero variance during fitting.", - UserWarning + f'Masked {num_zero_variance} data point(s) in reflectivity {i} due to zero variance during fitting.', + UserWarning, ) # Keep only points with non-zero variances diff --git a/src/easyreflectometry/model/model.py b/src/easyreflectometry/model/model.py index 6f7a203c..e4bdaac5 100644 --- a/src/easyreflectometry/model/model.py +++ b/src/easyreflectometry/model/model.py @@ -88,6 +88,7 @@ def __init__( scale = get_as_parameter('scale', scale, DEFAULTS) background = get_as_parameter('background', background, DEFAULTS) self.color = color + self._is_default = False super().__init__( name=name, @@ -138,6 +139,20 @@ def remove_assembly(self, index: int) -> None: if self.interface is not None: self.interface().remove_item_from_model(assembly_unique_name, self.unique_name) + @property + def is_default(self) -> bool: + """Whether this model was created as a default placeholder.""" + return self._is_default + + @is_default.setter + def is_default(self, value: bool) -> None: + """Set whether this model is a default placeholder. + + :param value: True if the model is a default placeholder. + :type value: bool + """ + self._is_default = value + @property def resolution_function(self) -> ResolutionFunction: """Return the resolution function.""" diff --git a/src/easyreflectometry/orso_utils.py b/src/easyreflectometry/orso_utils.py index 23756eba..494ed248 100644 --- a/src/easyreflectometry/orso_utils.py +++ b/src/easyreflectometry/orso_utils.py @@ -1,4 +1,5 @@ import logging +import warnings import numpy as np import scipp as sc @@ -19,23 +20,35 @@ logger = logging.getLogger(__name__) -def LoadOrso(orso_str: str): +def LoadOrso(orso_data): """Load a model from an ORSO file.""" - sample = load_orso_model(orso_str) - data = load_orso_data(orso_str) - + orso_obj = _coerce_orso_object(orso_data) + sample = load_orso_model(orso_obj) + data = load_orso_data(orso_obj) return sample, data + +def _coerce_orso_object(orso_input): + """Return a parsed ORSO object list from either a path or pre-parsed input.""" + try: + if orso_input and hasattr(orso_input[0], 'info'): + return orso_input + except (TypeError, IndexError): + pass + return orso.load_orso(orso_input) + + def load_data_from_orso_file(fname: str) -> sc.DataGroup: """Load data from an ORSO file.""" try: orso_data = orso.load_orso(fname) except Exception as e: - raise ValueError(f"Error loading ORSO file: {e}") + raise ValueError(f'Error loading ORSO file: {e}') return load_orso_data(orso_data) -def load_orso_model(orso_str: str) -> Sample: + +def load_orso_model(orso_data) -> Sample: """ Load a model from an ORSO file and return a Sample object. @@ -43,18 +56,29 @@ def load_orso_model(orso_str: str) -> Sample: as a simple "stack" string, e.g. 'air | m1 | SiO2 | Si'. This gets parsed by the ORSO library and converted into an ORSO Dataset object. - Args: - orso_str: The ORSO file content as a string - - Returns: - Sample: An EasyReflectometry Sample object + The stack is converted to a proper Sample structure: + - First layer -> Superphase assembly (thickness=0, roughness=0, both fixed) + - Middle layers -> 'Loaded layer' Multilayer assembly (parameters enabled) + - Last layer -> Subphase assembly (thickness=0 fixed, roughness enabled) - Raises: - ValueError: If ORSO layers could not be resolved + :param orso_data: Parsed ORSO dataset list (as returned by ``orso.load_orso``). + :type orso_data: list + :return: An EasyReflectometry Sample object. + :rtype: Sample + :raises ValueError: If ORSO layers could not be resolved or fewer than 2 layers. """ - # Extract stack string and create ORSO sample model - stack_str = orso_str[0].info.data_source.sample.model.stack - orso_sample = model_language.SampleModel(stack=stack_str) + # Extract stack string and layer definitions from ORSO sample model + sample_model = orso_data[0].info.data_source.sample.model + if sample_model is None: + warnings.warn( + 'ORSO file does not contain a sample model definition. Only experimental data can be loaded from this file.', + UserWarning, + stacklevel=2, + ) + return None + stack_str = sample_model.stack + layers_dict = sample_model.layers if hasattr(sample_model, 'layers') else None + orso_sample = model_language.SampleModel(stack=stack_str, layers=layers_dict) # Try to resolve layers using different methods try: @@ -64,9 +88,12 @@ def load_orso_model(orso_str: str) -> Sample: # Handle case where layers are not resolved correctly if not orso_layers: - raise ValueError("Could not resolve ORSO layers.") + raise ValueError('Could not resolve ORSO layers.') + + if len(orso_layers) < 2: + raise ValueError('ORSO stack must contain at least 2 layers (superphase and subphase).') - logger.debug(f"Resolved layers: {orso_layers}") + logger.debug(f'Resolved layers: {orso_layers}') # Convert ORSO layers to EasyReflectometry layers erl_layers = [] @@ -74,13 +101,34 @@ def load_orso_model(orso_str: str) -> Sample: erl_layer = _convert_orso_layer_to_erl(layer) erl_layers.append(erl_layer) - # Create a Multilayer object with the extracted layers - multilayer = Multilayer(erl_layers, name='Multi Layer Sample from ORSO') + # Create Superphase from first layer (thickness=0, roughness=0, both fixed) + superphase_layer = erl_layers[0] + superphase_layer.thickness.value = 0.0 + superphase_layer.roughness.value = 0.0 + superphase_layer.thickness.fixed = True + superphase_layer.roughness.fixed = True + superphase = Multilayer(superphase_layer, name='Superphase') + + # Create Subphase from last layer (thickness=0 fixed, roughness enabled) + subphase_layer = erl_layers[-1] + subphase_layer.thickness.value = 0.0 + subphase_layer.thickness.fixed = True + subphase_layer.roughness.fixed = False + subphase = Multilayer(subphase_layer, name='Subphase') # Create Sample from the file - sample_info = orso_str[0].info.data_source.sample + sample_info = orso_data[0].info.data_source.sample sample_name = sample_info.name if sample_info.name else 'ORSO Sample' - sample = Sample(multilayer, name=sample_name) + + # Build Sample based on number of layers + if len(erl_layers) == 2: + # Only superphase and subphase, no middle layers + sample = Sample(superphase, subphase, name=sample_name) + else: + # Create middle layer assembly from layers between first and last + middle_layers = erl_layers[1:-1] + loaded_layer = Multilayer(middle_layers, name='Loaded layer') + sample = Sample(superphase, loaded_layer, subphase, name=sample_name) return sample @@ -88,46 +136,75 @@ def load_orso_model(orso_str: str) -> Sample: def _convert_orso_layer_to_erl(layer): """Helper function to convert an ORSO layer to an EasyReflectometry layer""" material = layer.material - m_name = material.formula if material.formula is not None else layer.name + # Prefer original_name for material name, fall back to formula if available + m_name = layer.original_name if layer.original_name is not None else material.formula - # Get SLD values - m_sld, m_isld = _get_sld_values(material, m_name) + # Get SLD values (use formula for density calculation if available) + formula_for_calc = material.formula if material.formula is not None else m_name + m_sld, m_isld = _get_sld_values(material, formula_for_calc) # Create and return ERL layer return Layer( material=Material(sld=m_sld, isld=m_isld, name=m_name), thickness=layer.thickness.magnitude if layer.thickness is not None else 0.0, roughness=layer.roughness.magnitude if layer.roughness is not None else 0.0, - name=layer.original_name if layer.original_name is not None else m_name + name=layer.original_name if layer.original_name is not None else m_name, ) def _get_sld_values(material, material_name): - """Extract SLD values from material, calculating from density if needed""" + """Extract SLD values from material, calculating from density if needed + + Note: ORSO stores SLD in absolute units (A^-2), but the internal representation + uses 10^-6 A^-2. When reading directly from ORSO, we multiply by 1e6 to convert. + When calculating from mass density, MaterialDensity already returns the correct units. + """ if material.sld is None and material.mass_density is not None: # Calculate SLD from mass density + # MaterialDensity already returns values in 10^-6 A^-2 units m_density = material.mass_density.magnitude - density = MaterialDensity( - chemical_structure=material_name, - density=m_density - ) + density = MaterialDensity(chemical_structure=material_name, density=m_density) m_sld = density.sld.value m_isld = density.isld.value + elif material.sld is None: + # No SLD and no mass density available, default to 0.0 + m_sld = 0.0 + m_isld = 0.0 else: + # ORSO stores SLD in absolute units (A^-2) + # Convert to internal representation (10^-6 A^-2) by multiplying by 1e6 if isinstance(material.sld, ComplexValue): - m_sld = material.sld.real - m_isld = material.sld.imag + raw_sld = material.sld.real + m_sld = raw_sld * 1e6 + m_isld = material.sld.imag * 1e6 else: - m_sld = material.sld + raw_sld = material.sld + m_sld = raw_sld * 1e6 m_isld = 0.0 + if raw_sld != 0.0 and abs(raw_sld) > 1e-2: + warnings.warn( + f'ORSO SLD value {raw_sld} for "{material_name}" seems large for ' + f'absolute units (A^-2). Verify the file stores SLD in A^-2, not ' + f'10^-6 A^-2, as the value is multiplied by 1e6 internally.', + UserWarning, + stacklevel=3, + ) return m_sld, m_isld -def load_orso_data(orso_str: str) -> DataSet1D: + +def load_orso_data(orso_data) -> DataSet1D: + """Convert parsed ORSO dataset objects into a scipp DataGroup. + + :param orso_data: Parsed ORSO dataset list (as returned by ``orso.load_orso``). + :type orso_data: list + :return: A scipp DataGroup with data, coords, and attrs. + :rtype: sc.DataGroup + """ data = {} coords = {} attrs = {} - for i, o in enumerate(orso_str): + for i, o in enumerate(orso_data): name = i if o.info.data_set is not None: name = o.info.data_set diff --git a/src/easyreflectometry/project.py b/src/easyreflectometry/project.py index 196b9af1..13b4fdda 100644 --- a/src/easyreflectometry/project.py +++ b/src/easyreflectometry/project.py @@ -17,7 +17,10 @@ from easyreflectometry.calculators import CalculatorFactory from easyreflectometry.data import DataSet1D from easyreflectometry.data import load_as_dataset +from easyreflectometry.data.measurement import extract_orso_title +from easyreflectometry.data.measurement import load_data_from_orso_file from easyreflectometry.fitting import MultiFitter +from easyreflectometry.model import LinearSpline from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm @@ -70,11 +73,20 @@ def reset(self): @property def parameters(self) -> List[Parameter]: - """Get all parameters from all models in the project.""" + """Get all unique parameters from all models in the project. + + Parameters shared across multiple models (e.g. material SLD) are + only included once to avoid double-counting. + """ parameters = [] + seen_ids: set[int] = set() if self._models is not None: for model in self._models: - parameters.extend(model.get_parameters()) + for param in model.get_parameters(): + pid = id(param) + if pid not in seen_ids: + seen_ids.add(pid) + parameters.append(param) return parameters @property @@ -258,6 +270,8 @@ def load_orso_file(self, path: Union[Path, str]) -> None: model, data = LoadOrso(path) if model is not None: + if isinstance(model, Sample): + model = Model(sample=model, name=model.name) self.models = ModelCollection([model]) else: self.default_model() @@ -268,49 +282,148 @@ def load_orso_file(self, path: Union[Path, str]) -> None: self._with_experiments = True pass - def set_sample_from_orso(self, sample) -> None: + def set_sample_from_orso(self, sample: Sample) -> None: + """Replace the current project model collection with a single model built from an ORSO-parsed sample. + + This is a convenience helper for the ORSO import pipeline where a complete + :class:`~easyreflectometry.sample.Sample` is constructed elsewhere. + + :param sample: Sample to set as the project's (single) model. + :type sample: easyreflectometry.sample.Sample + :return: ``None``. + :rtype: None + """ model = Model(sample=sample) self.models = ModelCollection([model]) + def add_sample_from_orso(self, sample: Sample) -> None: + """Add a new model with the given sample to the existing model collection. + + The created model is appended to :attr:`models`, its calculator interface is + set to the project's current calculator, and any materials referenced in the + sample are added to the project's material collection. + + After adding the model, :attr:`current_model_index` is updated to point to + the newly added model. + + :param sample: Sample to add as a new model. + :type sample: easyreflectometry.sample.Sample + :return: ``None``. + :rtype: None + """ + if sample is None: + raise ValueError('The ORSO file does not contain a valid sample model definition.') + model = Model(sample=sample) + self.models.add_model(model) + # Set interface after adding to collection + model.interface = self._calculator + # Extract materials from the new model and add to project materials + self._materials.extend(self._get_materials_from_model(model)) + # Switch to the newly added model so its data is visible in the UI + self.current_model_index = len(self._models) - 1 + + def replace_models_from_orso(self, sample: Sample) -> None: + """Replace all models and materials with a single model from an ORSO sample. + + All existing models and their associated materials are removed. A new + model is created from *sample*, assigned to the project's calculator, + and the material collection is rebuilt from the new model only. + + :param sample: Sample to set as the project's only model. + :type sample: easyreflectometry.sample.Sample + :return: ``None``. + :rtype: None + """ + if sample is None: + raise ValueError('The ORSO file does not contain a valid sample model definition.') + model = Model(sample=sample) + if sample.name: + model.user_data['original_name'] = sample.name # Store original name for reference + self.models = ModelCollection([model]) + model.interface = self._calculator + self._materials = self._get_materials_from_model(model) + self.current_model_index = 0 + + def _get_materials_from_model(self, model: Model) -> 'MaterialCollection': + """Get all materials from a single model's sample.""" + materials_in_model = MaterialCollection(populate_if_none=False) + for assembly in model.sample: + for layer in assembly.layers: + if layer.material not in materials_in_model: + materials_in_model.append(layer.material) + return materials_in_model + + def _apply_experiment_metadata( + self, + path: Union[Path, str], + experiment: DataSet1D, + fallback_name: str, + ) -> None: + """Set experiment name from ORSO title and configure the resolution function. + + :param path: Path to the experiment data file. + :param experiment: The loaded experiment dataset to configure. + :param fallback_name: Name to use when no ORSO title is available. + """ + # Prefer ORSO title when available (keeps UI descriptive) + title = None + try: + data_group = load_data_from_orso_file(str(path)) + data_key = list(data_group['data'].keys())[0] + title = extract_orso_title(data_group, data_key) + except (KeyError, AttributeError, ValueError, IndexError): + title = None + + if title: + experiment.name = title + elif not experiment.name or experiment.name == 'Series': + experiment.name = fallback_name + + def _apply_resolution_function( + self, + experiment: DataSet1D, + model: Model, + ) -> None: + """Set the resolution function on *model* based on variance data in *experiment*. + + Prefers Pointwise when q-resolution (xe) data is present, otherwise falls + back to LinearSpline when reflectivity error (ye) data is present. + + :param experiment: The experiment whose variance data drives the choice. + :param model: The model whose resolution function is set. + """ + if sum(experiment.xe) != 0: + resolution_function = Pointwise(q_data_points=[experiment.x, experiment.y, experiment.xe]) + model.resolution_function = resolution_function + elif sum(experiment.ye) != 0: + resolution_function = LinearSpline( + q_data_points=experiment.x, + fwhm_values=np.sqrt(experiment.ye), + ) + model.resolution_function = resolution_function + def load_new_experiment(self, path: Union[Path, str]) -> None: new_experiment = load_as_dataset(str(path)) new_index = len(self._experiments) - new_experiment.name = f'Experiment {new_index}' + model_index = 0 if new_index < len(self.models): model_index = new_index + + self._apply_experiment_metadata(path, new_experiment, f'Experiment {new_index}') new_experiment.model = self.models[model_index] self._experiments[new_index] = new_experiment - # self._current_model_index = new_index self._with_experiments = True - - # Set the resolution function if variance data is present - if sum(new_experiment.ye) != 0: - q = new_experiment.x - reflectivity = new_experiment.y - q_error = new_experiment.xe - # TODO: set resolution function based on value of control in GUI - resolution_function = Pointwise(q_data_points=[q, reflectivity, q_error]) - self.models[model_index].resolution_function = resolution_function + self._apply_resolution_function(new_experiment, self.models[model_index]) def load_experiment_for_model_at_index(self, path: Union[Path, str], index: Optional[int] = 0) -> None: - self._experiments[index] = load_as_dataset(str(path)) - self._experiments[index].name = f'Experiment {index}' - self._experiments[index].model = self.models[index] + experiment = load_as_dataset(str(path)) + self._apply_experiment_metadata(path, experiment, f'Experiment {index}') + experiment.model = self.models[index] + self._experiments[index] = experiment self._with_experiments = True - - # Set the resolution function if variance data is present - if sum(self._experiments[index].ye) != 0: - q = self._experiments[index].x - reflectivity = self._experiments[index].y - q_error = self._experiments[index].xe - resolution_function = Pointwise(q_data_points=[q, reflectivity, q_error]) - # resolution_function = LinearSpline( - # q_data_points=self._experiments[index].y, - # fwhm_values=np.sqrt(self._experiments[index].ye), - # ) - self._models[index].resolution_function = resolution_function + self._apply_resolution_function(experiment, self._models[index]) def sld_data_for_model_at_index(self, index: int = 0) -> DataSet1D: self.models[index].interface = self._calculator @@ -361,8 +474,68 @@ def default_model(self): ] sample = Sample(*assemblies, interface=self._calculator) model = Model(sample=sample, interface=self._calculator) + model.is_default = True self.models = ModelCollection([model]) + def is_default_model(self, index: int) -> bool: + """Check if the model at the given index is a default model. + + :param index: Index of the model to check. + :type index: int + :return: True if the model was created as a default placeholder. + :rtype: bool + """ + if index < 0 or index >= len(self._models): + return False + + return self._models[index].is_default + + def remove_model_at_index(self, index: int) -> None: + """Remove the model at the given index. + + Removes the model from the model collection, removes the experiment at the + same index (if any), and reindexes experiments above the removed index so + model/experiment indices stay aligned. + + Adjusts the current model index if necessary. + + :param index: Index of the model to remove. + :type index: int + :raises IndexError: If the index is out of range. + :raises ValueError: If trying to remove the last remaining model. + """ + if index < 0 or index >= len(self._models): + raise IndexError(f'Model index {index} out of range') + + if len(self._models) <= 1: + raise ValueError('Cannot remove the last model from the project') + + # Remove the model from the collection + self._models.pop(index) + + # Remove experiment mapped to the removed model index. + if index in self._experiments: + self._experiments.pop(index) + + # Reindex experiments above the removed model index to keep mapping aligned. + reindexed_experiments: dict[int, DataSet1D] = {} + for exp_index, experiment in sorted(self._experiments.items()): + if exp_index > index: + reindexed_experiments[exp_index - 1] = experiment + else: + reindexed_experiments[exp_index] = experiment + self._experiments = reindexed_experiments + + # Adjust current model index if necessary + if self._current_model_index >= len(self._models): + self._current_model_index = len(self._models) - 1 + elif self._current_model_index > index: + self._current_model_index -= 1 + + # Reset assembly and layer indices for the new current model + self._current_assembly_index = 0 + self._current_layer_index = 0 + def add_material(self, material: MaterialCollection) -> None: if material in self._materials: print(f'WARNING: Material {material} is already in material collection') diff --git a/tests/calculators/refl1d/test_refl1d_calculator.py b/tests/calculators/refl1d/test_refl1d_calculator.py index 4dff8a9b..ba8c8d35 100644 --- a/tests/calculators/refl1d/test_refl1d_calculator.py +++ b/tests/calculators/refl1d/test_refl1d_calculator.py @@ -61,7 +61,7 @@ def test_reflectity_profile(self): 5.7605e-07, 2.3775e-07, 1.3093e-07, - 1.0520e-07 + 1.0520e-07, ] assert_almost_equal(p.reflectity_profile(q, 'MyModel'), expected, decimal=4) @@ -106,7 +106,7 @@ def test_calculate2(self): 1.0968e-06, 4.5635e-07, 3.4120e-07, - 2.7505e-07 + 2.7505e-07, ] assert_almost_equal(actual, expected, decimal=4) diff --git a/tests/calculators/refl1d/test_refl1d_wrapper.py b/tests/calculators/refl1d/test_refl1d_wrapper.py index 725aca6a..e19dfe42 100644 --- a/tests/calculators/refl1d/test_refl1d_wrapper.py +++ b/tests/calculators/refl1d/test_refl1d_wrapper.py @@ -232,7 +232,7 @@ def test_calculate(self): 5.7605e-07, 2.3775e-07, 1.3093e-07, - 1.0520e-07 + 1.0520e-07, ] assert_almost_equal(p.calculate(q, 'MyModel'), expected, decimal=4) @@ -276,7 +276,7 @@ def test_calculate_three_items(self): 1.0968e-06, 4.5635e-07, 3.4120e-07, - 2.7505e-07 + 2.7505e-07, ] assert_almost_equal(p.calculate(q, 'MyModel'), expected, decimal=4) @@ -396,7 +396,7 @@ def test_get_polarized_probe_oversampling(): probe = _get_polarized_probe(q_array=q, dq_array=dq, model_name=model_name, storage=storage, oversampling_factor=2) # Then - assert len(probe.xs[0].calc_Qo) == 2*len(q) + assert len(probe.xs[0].calc_Qo) == 2 * len(q) def test_get_polarized_probe_polarization(): diff --git a/tests/data/test_data_store.py b/tests/data/test_data_store.py index 17f837e2..84acf4f8 100644 --- a/tests/data/test_data_store.py +++ b/tests/data/test_data_store.py @@ -29,13 +29,7 @@ def test_constructor_default_values(self): def test_constructor_with_values(self): # When data = DataSet1D( - x=[1, 2, 3], - y=[4, 5, 6], - ye=[7, 8, 9], - xe=[10, 11, 12], - x_label='label_x', - y_label='label_y', - name='MyDataSet1D' + x=[1, 2, 3], y=[4, 5, 6], ye=[7, 8, 9], xe=[10, 11, 12], x_label='label_x', y_label='label_y', name='MyDataSet1D' ) # Then @@ -116,9 +110,7 @@ def test_is_simulation_property(self): def test_data_points(self): # When - data = DataSet1D( - x=[1, 2, 3], y=[4, 5, 6], ye=[7, 8, 9], xe=[10, 11, 12] - ) + data = DataSet1D(x=[1, 2, 3], y=[4, 5, 6], ye=[7, 8, 9], xe=[10, 11, 12]) # Then points = list(data.data_points()) @@ -126,9 +118,7 @@ def test_data_points(self): def test_repr(self): # When - data = DataSet1D( - x=[1, 2, 3], y=[4, 5, 6], x_label='Q', y_label='R' - ) + data = DataSet1D(x=[1, 2, 3], y=[4, 5, 6], x_label='Q', y_label='R') # Then expected = "1D DataStore of 'Q' Vs 'R' with 3 data points" @@ -194,7 +184,7 @@ def test_setitem(self): item1 = DataSet1D(name='item1') item2 = DataSet1D(name='item2') store = DataStore(item1) - + # When store[0] = item2 @@ -314,4 +304,3 @@ def test_constructor_with_custom_datastores(self): assert project.sim_data == sim_store assert project.exp_data.name == 'CustomExp' assert project.sim_data.name == 'CustomSim' - diff --git a/tests/summary/test_summary.py b/tests/summary/test_summary.py index 319a1c9b..65de9ba3 100644 --- a/tests/summary/test_summary.py +++ b/tests/summary/test_summary.py @@ -129,7 +129,7 @@ def test_experiments_section(self, project: Project) -> None: html = summary._experiments_section() # Expect - assert 'Experiment 0' in html + assert 'Example data file from refnx docs' in html assert 'No. of data points' in html assert '408' in html assert 'Resolution function' in html @@ -177,7 +177,7 @@ def test_save_sld_plot(self, project: Project, tmp_path) -> None: # Expect assert os.path.exists(file_path) - @pytest.mark.skip(reason="Matplotlib issue with headless CI environments") + @pytest.mark.skip(reason='Matplotlib issue with headless CI environments') def test_save_fit_experiment_plot(self, project: Project, tmp_path) -> None: # When summary = Summary(project) diff --git a/tests/test_data.py b/tests/test_data.py index c16aa259..0ee95d94 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -224,8 +224,7 @@ def test_load_txt_three_columns(self): assert coords_name in er_data['coords'] # xe should be zeros for 3-column file - assert_almost_equal(er_data['coords'][coords_name].variances, - np.zeros_like(er_data['coords'][coords_name].values)) + assert_almost_equal(er_data['coords'][coords_name].variances, np.zeros_like(er_data['coords'][coords_name].values)) def test_load_txt_with_zero_errors(self): fpath = os.path.join(PATH_STATIC, 'ref_zero_var.txt') @@ -246,6 +245,7 @@ def test_load_txt_file_not_found(self): def test_load_txt_insufficient_columns(self): # Create a temporary file with insufficient columns import tempfile + with tempfile.NamedTemporaryFile(mode='w', suffix='.txt', delete=False) as f: f.write('1.0 2.0\n') # Only 2 columns temp_path = f.name @@ -272,7 +272,7 @@ def test_load_orso_multiple_datasets(self): if data_key.replace('R_', '') in coord_key: coord_key_found = True break - assert coord_key_found, f"No corresponding coord found for {data_key}" + assert coord_key_found, f'No corresponding coord found for {data_key}' def test_load_orso_with_attrs(self): fpath = os.path.join(PATH_STATIC, 'test_example1.ort') diff --git a/tests/test_fitting.py b/tests/test_fitting.py index 0b02ed82..fdeba385 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -86,7 +86,7 @@ def test_fitting_with_zero_variance(): # First, load the raw data to count zero variance points raw_data = np.loadtxt(fpath, delimiter=',', comments='#') zero_variance_count = np.sum(raw_data[:, 2] == 0.0) # Error column - assert zero_variance_count == 6, f"Expected 6 zero variance points, got {zero_variance_count}" + assert zero_variance_count == 6, f'Expected 6 zero variance points, got {zero_variance_count}' # Load data through the measurement module (which already filters zero variance) data = load(fpath) @@ -129,12 +129,11 @@ def test_fitting_with_zero_variance(): # Capture warnings during fitting - check if zero variance points still exist in the data # and are properly handled by the fitting method with warnings.catch_warnings(record=True) as w: - warnings.simplefilter("always") + warnings.simplefilter('always') analysed = fitter.fit(data) # Check if any zero variance warnings were issued during fitting - fitting_warnings = [str(warning.message) for warning in w - if "zero variance during fitting" in str(warning.message)] + fitting_warnings = [str(warning.message) for warning in w if 'zero variance during fitting' in str(warning.message)] # The fitting method should handle zero variance points gracefully # If there are any zero variance points remaining in the data, they should be masked @@ -142,15 +141,15 @@ def test_fitting_with_zero_variance(): if len(fitting_warnings) > 0: # Verify the warning message format and that it mentions masking points for warning_msg in fitting_warnings: - assert "Masked" in warning_msg and "zero variance during fitting" in warning_msg - print(f"Info: {warning_msg}") # Log for debugging + assert 'Masked' in warning_msg and 'zero variance during fitting' in warning_msg + print(f'Info: {warning_msg}') # Log for debugging # Basic checks that fitting completed # The keys will be based on the filename, not just '0' model_keys = [k for k in analysed.keys() if k.endswith('_model')] sld_keys = [k for k in analysed.keys() if k.startswith('SLD_')] - assert len(model_keys) > 0, f"No model keys found in {list(analysed.keys())}" - assert len(sld_keys) > 0, f"No SLD keys found in {list(analysed.keys())}" + assert len(model_keys) > 0, f'No model keys found in {list(analysed.keys())}' + assert len(sld_keys) > 0, f'No SLD keys found in {list(analysed.keys())}' assert 'success' in analysed.keys() @@ -172,14 +171,12 @@ def test_fitting_with_manual_zero_variance(): variances[30:32] = 0.0 # 2 more zero variance points # Create scipp DataGroup manually - data = sc.DataGroup({ - 'coords': { - 'Qz_0': sc.array(dims=['Qz_0'], values=qz_values) - }, - 'data': { - 'R_0': sc.array(dims=['Qz_0'], values=r_values, variances=variances) + data = sc.DataGroup( + { + 'coords': {'Qz_0': sc.array(dims=['Qz_0'], values=qz_values)}, + 'data': {'R_0': sc.array(dims=['Qz_0'], values=r_values, variances=variances)}, } - }) + ) # Create a simple model for fitting si = Material(2.07, 0, 'Si') @@ -214,16 +211,15 @@ def test_fitting_with_manual_zero_variance(): # Capture warnings during fitting with warnings.catch_warnings(record=True) as w: - warnings.simplefilter("always") + warnings.simplefilter('always') analysed = fitter.fit(data) # Check that warnings were issued about zero variance points - fitting_warnings = [str(warning.message) for warning in w - if "zero variance during fitting" in str(warning.message)] + fitting_warnings = [str(warning.message) for warning in w if 'zero variance during fitting' in str(warning.message)] # Should have one warning about the 7 zero variance points (5 + 2) - assert len(fitting_warnings) == 1, f"Expected 1 warning, got {len(fitting_warnings)}: {fitting_warnings}" - assert "Masked 7 data point(s)" in fitting_warnings[0], f"Unexpected warning content: {fitting_warnings[0]}" + assert len(fitting_warnings) == 1, f'Expected 1 warning, got {len(fitting_warnings)}: {fitting_warnings}' + assert 'Masked 7 data point(s)' in fitting_warnings[0], f'Unexpected warning content: {fitting_warnings[0]}' # Basic checks that fitting completed despite zero variance points assert 'R_0_model' in analysed.keys() assert 'SLD_0' in analysed.keys() diff --git a/tests/test_measurement_comprehensive.py b/tests/test_measurement_comprehensive.py index 3842fbb5..e9bf6ffe 100644 --- a/tests/test_measurement_comprehensive.py +++ b/tests/test_measurement_comprehensive.py @@ -33,7 +33,7 @@ def test_load_function_with_orso_file(self): """Test that load() correctly identifies and loads ORSO files.""" fpath = os.path.join(PATH_STATIC, 'test_example1.ort') result = load(fpath) - + assert 'data' in result assert 'coords' in result assert len(result['data']) > 0 @@ -43,7 +43,7 @@ def test_load_function_with_txt_file(self): """Test that load() falls back to txt loading for non-ORSO files.""" fpath = os.path.join(PATH_STATIC, 'test_example1.txt') result = load(fpath) - + assert 'data' in result assert 'coords' in result assert 'R_test_example1' in result['data'] @@ -53,7 +53,7 @@ def test_load_as_dataset_returns_dataset1d(self): """Test that load_as_dataset returns a proper DataSet1D object.""" fpath = os.path.join(PATH_STATIC, 'test_example1.txt') dataset = load_as_dataset(fpath) - + assert isinstance(dataset, DataSet1D) assert hasattr(dataset, 'x') assert hasattr(dataset, 'y') @@ -65,7 +65,7 @@ def test_load_as_dataset_extracts_correct_basename(self): """Test that load_as_dataset correctly extracts file basename.""" fpath = os.path.join(PATH_STATIC, 'ref_concat_1.txt') dataset = load_as_dataset(fpath) - + # Should work without error and have data assert len(dataset.x) > 0 assert len(dataset.y) > 0 @@ -74,12 +74,12 @@ def test_merge_datagroups_preserves_all_data(self): """Test that merge_datagroups combines multiple data groups correctly.""" fpath1 = os.path.join(PATH_STATIC, 'test_example1.txt') fpath2 = os.path.join(PATH_STATIC, 'ref_concat_1.txt') - + group1 = load(fpath1) group2 = load(fpath2) - + merged = merge_datagroups(group1, group2) - + # Should have data from both groups assert len(merged['data']) >= len(group1['data']) assert len(merged['coords']) >= len(group1['coords']) @@ -88,9 +88,9 @@ def test_merge_datagroups_single_group(self): """Test that merge_datagroups works with a single group.""" fpath = os.path.join(PATH_STATIC, 'test_example1.ort') group = load(fpath) - + merged = merge_datagroups(group) - + # Should be equivalent to original assert len(merged['data']) == len(group['data']) assert len(merged['coords']) == len(group['coords']) @@ -99,7 +99,7 @@ def test_load_txt_handles_comma_delimiter(self): """Test that _load_txt correctly handles comma-delimited files.""" fpath = os.path.join(PATH_STATIC, 'ref_concat_1.txt') result = _load_txt(fpath) - + assert 'data' in result assert 'coords' in result # Should successfully parse comma-delimited data @@ -110,11 +110,10 @@ def test_load_txt_handles_three_columns(self): """Test that _load_txt handles files with only 3 columns (no xe).""" fpath = os.path.join(PATH_STATIC, 'ref_concat_1.txt') result = _load_txt(fpath) - + coords_key = list(result['coords'].keys())[0] # xe should be zeros - assert_array_equal(result['coords'][coords_key].variances, - np.zeros_like(result['coords'][coords_key].values)) + assert_array_equal(result['coords'][coords_key].variances, np.zeros_like(result['coords'][coords_key].values)) def test_load_txt_with_insufficient_columns(self): """Test that _load_txt raises error for files with too few columns.""" @@ -122,7 +121,7 @@ def test_load_txt_with_insufficient_columns(self): with tempfile.NamedTemporaryFile(mode='w', suffix='.txt', delete=False) as f: f.write('1.0 2.0\n3.0 4.0\n') temp_path = f.name - + try: with pytest.raises(ValueError, match='File must contain at least 3 columns'): _load_txt(temp_path) @@ -133,7 +132,7 @@ def test_load_orso_with_multiple_datasets(self): """Test that _load_orso handles files with multiple datasets.""" fpath = os.path.join(PATH_STATIC, 'test_example2.ort') result = load_data_from_orso_file(fpath) - + # Should have multiple data entries assert len(result['data']) > 1 assert 'attrs' in result @@ -142,7 +141,7 @@ def test_load_orso_preserves_metadata(self): """Test that _load_orso preserves ORSO metadata in attrs.""" fpath = os.path.join(PATH_STATIC, 'test_example1.ort') result = load_data_from_orso_file(fpath) - + assert 'attrs' in result # Should have orso_header in attrs for data_key in result['data']: @@ -159,15 +158,9 @@ def test_constructor_all_parameters(self): y = [10, 20, 30, 40] xe = [0.1, 0.1, 0.1, 0.1] ye = [1, 2, 3, 4] - - dataset = DataSet1D( - name='TestData', - x=x, y=y, xe=xe, ye=ye, - x_label='Q (Å⁻¹)', - y_label='Reflectivity', - model=None - ) - + + dataset = DataSet1D(name='TestData', x=x, y=y, xe=xe, ye=ye, x_label='Q (Å⁻¹)', y_label='Reflectivity', model=None) + assert dataset.name == 'TestData' assert_array_equal(dataset.x, np.array(x)) assert_array_equal(dataset.y, np.array(y)) @@ -182,7 +175,7 @@ def test_is_experiment_vs_simulation_properties(self): sim_data = DataSet1D(x=[1, 2], y=[3, 4]) assert sim_data.is_simulation is True assert sim_data.is_experiment is False - + # Dataset with model is experiment exp_data = DataSet1D(x=[1, 2], y=[3, 4], model=Mock()) assert exp_data.is_experiment is True @@ -190,13 +183,8 @@ def test_is_experiment_vs_simulation_properties(self): def test_data_points_iterator(self): """Test the data_points method returns correct tuples.""" - dataset = DataSet1D( - x=[1, 2, 3], - y=[10, 20, 30], - xe=[0.1, 0.2, 0.3], - ye=[1, 2, 3] - ) - + dataset = DataSet1D(x=[1, 2, 3], y=[10, 20, 30], xe=[0.1, 0.2, 0.3], ye=[1, 2, 3]) + points = list(dataset.data_points()) expected = [(1, 10, 1, 0.1), (2, 20, 2, 0.2), (3, 30, 3, 0.3)] assert points == expected @@ -205,20 +193,15 @@ def test_model_property_with_background_setting(self): """Test that setting model updates background to minimum y value.""" dataset = DataSet1D(x=[1, 2, 3, 4], y=[5, 1, 8, 3]) mock_model = Mock() - + dataset.model = mock_model - + assert mock_model.background == 1 # minimum of [5, 1, 8, 3] def test_repr_string_representation(self): """Test the string representation of DataSet1D.""" - dataset = DataSet1D( - x=[1, 2, 3], - y=[4, 5, 6], - x_label='Momentum Transfer', - y_label='Intensity' - ) - + dataset = DataSet1D(x=[1, 2, 3], y=[4, 5, 6], x_label='Momentum Transfer', y_label='Intensity') + expected = "1D DataStore of 'Momentum Transfer' Vs 'Intensity' with 3 data points" assert str(dataset) == expected @@ -230,19 +213,19 @@ def test_datastore_as_sequence(self): """Test DataStore behaves like a sequence.""" item1 = DataSet1D(name='item1', x=[1], y=[2]) item2 = DataSet1D(name='item2', x=[3], y=[4]) - + store = DataStore(item1, item2, name='TestStore') - + # Test sequence operations assert len(store) == 2 assert store[0].name == 'item1' assert store[1].name == 'item2' - + # Test item replacement item3 = DataSet1D(name='item3', x=[5], y=[6]) store[0] = item3 assert store[0].name == 'item3' - + # Test deletion del store[0] assert len(store) == 1 @@ -254,12 +237,12 @@ def test_datastore_experiments_and_simulations_filtering(self): exp2 = DataSet1D(name='exp2', x=[3], y=[4], model=Mock()) sim1 = DataSet1D(name='sim1', x=[5], y=[6]) sim2 = DataSet1D(name='sim2', x=[7], y=[8]) - + store = DataStore(exp1, sim1, exp2, sim2) - + experiments = store.experiments simulations = store.simulations - + assert len(experiments) == 2 assert len(simulations) == 2 assert all(item.is_experiment for item in experiments) @@ -269,9 +252,9 @@ def test_datastore_append_method(self): """Test append method adds items correctly.""" store = DataStore() item = DataSet1D(name='new_item', x=[1], y=[2]) - + store.append(item) - + assert len(store) == 1 assert store[0] == item @@ -282,7 +265,7 @@ class TestProjectDataComprehensive: def test_project_data_initialization(self): """Test ProjectData initializes with correct default values.""" project = ProjectData() - + assert project.name == 'DataStore' assert isinstance(project.exp_data, DataStore) assert isinstance(project.sim_data, DataStore) @@ -293,13 +276,9 @@ def test_project_data_with_custom_stores(self): """Test ProjectData with custom experiment and simulation stores.""" custom_exp = DataStore(name='CustomExp') custom_sim = DataStore(name='CustomSim') - - project = ProjectData( - name='MyProject', - exp_data=custom_exp, - sim_data=custom_sim - ) - + + project = ProjectData(name='MyProject', exp_data=custom_exp, sim_data=custom_sim) + assert project.name == 'MyProject' assert project.exp_data == custom_exp assert project.sim_data == custom_sim @@ -307,13 +286,13 @@ def test_project_data_with_custom_stores(self): def test_project_data_stores_independence(self): """Test that exp_data and sim_data are independent stores.""" project = ProjectData() - + exp_item = DataSet1D(name='exp', x=[1], y=[2], model=Mock()) sim_item = DataSet1D(name='sim', x=[3], y=[4]) - + project.exp_data.append(exp_item) project.sim_data.append(sim_item) - + assert len(project.exp_data) == 1 assert len(project.sim_data) == 1 assert project.exp_data[0] != project.sim_data[0] @@ -327,11 +306,11 @@ def test_complete_workflow_orso_file(self): # Load file fpath = os.path.join(PATH_STATIC, 'test_example1.ort') dataset = load_as_dataset(fpath) - + # Create project and add to experimental data project = ProjectData(name='MyAnalysis') project.exp_data.append(dataset) - + # Verify workflow assert len(project.exp_data) == 1 assert project.exp_data[0] == dataset @@ -342,11 +321,11 @@ def test_complete_workflow_txt_file(self): # Load file fpath = os.path.join(PATH_STATIC, 'ref_concat_1.txt') dataset = load_as_dataset(fpath) - + # Create project and add to simulation data (no model) project = ProjectData(name='MySimulation') project.sim_data.append(dataset) - + # Verify workflow assert len(project.sim_data) == 1 assert project.sim_data[0] == dataset @@ -357,13 +336,13 @@ def test_merge_multiple_files_workflow(self): # Load multiple files fpath1 = os.path.join(PATH_STATIC, 'test_example1.txt') fpath2 = os.path.join(PATH_STATIC, 'ref_concat_1.txt') - + group1 = load(fpath1) group2 = load(fpath2) - + # Merge data groups merged = merge_datagroups(group1, group2) - + # Create datasets from merged data # This tests that merged data can be used to create datasets assert len(merged['data']) >= 2 @@ -374,13 +353,13 @@ def test_error_handling_robustness(self): # Test mismatched array lengths with pytest.raises(ValueError, match='x and y must be the same length'): DataSet1D(x=[1, 2, 3], y=[4, 5]) - + # Test empty DataStore operations empty_store = DataStore() assert len(empty_store) == 0 assert len(empty_store.experiments) == 0 assert len(empty_store.simulations) == 0 - + # Test file not found with pytest.raises(FileNotFoundError): _load_txt('nonexistent_file.txt') @@ -391,14 +370,14 @@ def test_data_consistency_checks(self): original_x = [1, 2, 3, 4] original_y = [10, 20, 30, 40] dataset = DataSet1D(x=original_x, y=original_y) - + # Store in datastore store = DataStore(dataset) - + # Add to project project = ProjectData() project.sim_data = store - + # Verify data consistency retrieved_dataset = project.sim_data[0] assert_array_equal(retrieved_dataset.x, np.array(original_x)) @@ -407,4 +386,4 @@ def test_data_consistency_checks(self): if __name__ == '__main__': # Run all tests if script is executed directly - pytest.main([__file__, '-v']) \ No newline at end of file + pytest.main([__file__, '-v']) diff --git a/tests/test_orso_utils.py b/tests/test_orso_utils.py index 89dd07db..ebb662a1 100644 --- a/tests/test_orso_utils.py +++ b/tests/test_orso_utils.py @@ -2,12 +2,15 @@ # Copyright (c) 2025 DMSC import os +import warnings +from types import SimpleNamespace import pytest from orsopy.fileio import orso import easyreflectometry from easyreflectometry.orso_utils import LoadOrso +from easyreflectometry.orso_utils import _get_sld_values from easyreflectometry.orso_utils import load_data_from_orso_file from easyreflectometry.orso_utils import load_orso_data from easyreflectometry.orso_utils import load_orso_model @@ -18,14 +21,47 @@ @pytest.fixture def orso_data(): """Load the test ORSO data from Ni_example.ort.""" - return orso.load_orso(os.path.join(PATH_STATIC, "Ni_example.ort")) + return orso.load_orso(os.path.join(PATH_STATIC, 'Ni_example.ort')) def test_load_orso_model(orso_data): """Test loading a model from ORSO data.""" sample = load_orso_model(orso_data) assert sample is not None - assert sample.name == "Ni on Si" # Based on the file + assert sample.name == 'Ni on Si' # Based on the file + + # Verify sample structure: Superphase, Loaded layer, Subphase + # Stack in file: air | m1 | SiO2 | Si + assert len(sample) == 3 + + # Check Superphase (first layer from stack: air) + superphase = sample[0] + assert superphase.name == 'Superphase' + assert len(superphase.layers) == 1 + assert superphase.layers[0].material.name == 'air' + assert superphase.layers[0].thickness.value == 0.0 + assert superphase.layers[0].roughness.value == 0.0 + assert superphase.layers[0].thickness.fixed is True + assert superphase.layers[0].roughness.fixed is True + + # Check Loaded layer (middle layers: m1, SiO2) + loaded_layer = sample[1] + assert loaded_layer.name == 'Loaded layer' + assert len(loaded_layer.layers) == 2 + assert loaded_layer.layers[0].material.name == 'm1' # Uses original_name, not formula + assert loaded_layer.layers[0].thickness.value == 1000.0 # From layer definition + assert loaded_layer.layers[1].material.name == 'SiO2' + assert loaded_layer.layers[1].thickness.value == 10.0 # From layer definition + + # Check Subphase (last layer from stack: Si) + subphase = sample[2] + assert subphase.name == 'Subphase' + assert len(subphase.layers) == 1 + assert subphase.layers[0].material.name == 'Si' + assert subphase.layers[0].thickness.value == 0.0 + assert subphase.layers[0].thickness.fixed is True + # Subphase roughness should be enabled (not fixed) + assert subphase.layers[0].roughness.fixed is False def test_load_orso_data(orso_data): @@ -33,7 +69,7 @@ def test_load_orso_data(orso_data): data = load_orso_data(orso_data) assert data is not None # Check structure, e.g., has R_0 in data - assert "R_0" in data["data"] + assert 'R_0' in data['data'] def test_LoadOrso(orso_data): @@ -46,8 +82,93 @@ def test_LoadOrso(orso_data): def test_load_data_from_orso_file(): """Test loading data from ORSO file.""" - data = load_data_from_orso_file(os.path.join(PATH_STATIC, "Ni_example.ort")) + data = load_data_from_orso_file(os.path.join(PATH_STATIC, 'Ni_example.ort')) assert data is not None # Check it's a sc.DataGroup import scipp as sc - assert isinstance(data, sc.DataGroup) \ No newline at end of file + + assert isinstance(data, sc.DataGroup) + + +def test_orso_sld_unit_conversion(orso_data): + """Test that SLD values from ORSO are correctly converted from A^-2 to 10^-6 A^-2. + + ORSO stores SLD in absolute units (A^-2), e.g., 3.47e-06. + The internal representation uses 10^-6 A^-2, so the value should be 3.47. + """ + sample = load_orso_model(orso_data) + + # Check SiO2 layer (second layer in Loaded layer assembly) + # ORSO file has: sld: {real: 3.4700000000000002e-06, imag: 0.0} + # Expected internal value: 3.47 + loaded_layer = sample[1] + sio2_layer = loaded_layer.layers[1] + assert sio2_layer.material.name == 'SiO2' + assert abs(sio2_layer.material.sld.value - 3.47) < 1e-6, ( + f'Expected SLD ~3.47 (10^-6 A^-2), got {sio2_layer.material.sld.value}' + ) + + # Check Si subphase layer + # ORSO file has: sld: {real: 2.0699999999999997e-06, imag: 0.0} + # Expected internal value: 2.07 + subphase = sample[2] + si_layer = subphase.layers[0] + assert si_layer.material.name == 'Si' + assert abs(si_layer.material.sld.value - 2.07) < 1e-6, ( + f'Expected SLD ~2.07 (10^-6 A^-2), got {si_layer.material.sld.value}' + ) + + # Check air superphase layer + # ORSO file has: sld: {real: 0.0, imag: 0.0} + # Expected internal value: 0.0 + superphase = sample[0] + air_layer = superphase.layers[0] + assert air_layer.material.name == 'air' + assert abs(air_layer.material.sld.value - 0.0) < 1e-6, f'Expected SLD 0.0 (10^-6 A^-2), got {air_layer.material.sld.value}' + + +def test_LoadOrso_returns_two_items(orso_data): + """LoadOrso should return exactly two values: (sample, data).""" + result = LoadOrso(orso_data) + assert isinstance(result, tuple) + assert len(result) == 2 + sample, data = result + assert sample is not None + assert data is not None + + +def test_LoadOrso_with_invalid_file(tmp_path): + """LoadOrso should raise for a corrupt / non-ORSO file.""" + bad_file = tmp_path / 'bad.ort' + bad_file.write_text('this is not valid ORSO data') + with pytest.raises((ValueError, Exception)): + LoadOrso(str(bad_file)) + + +def test_LoadOrso_with_nonexistent_file(): + """LoadOrso should raise for a path that does not exist.""" + with pytest.raises((FileNotFoundError, ValueError, Exception)): + LoadOrso('/nonexistent/path/to/file.ort') + + +def test_get_sld_values_defaults_to_zero_when_sld_and_density_missing(): + """_get_sld_values should return (0.0, 0.0) when both sld and mass_density are None.""" + material = SimpleNamespace(sld=None, mass_density=None) + m_sld, m_isld = _get_sld_values(material, 'Unknown') + assert m_sld == 0.0 + assert m_isld == 0.0 + + +def test_load_orso_model_returns_none_and_warns_when_no_sample_model(): + """load_orso_model should return None and emit a warning when the ORSO file has no sample model.""" + orso_data = orso.load_orso(os.path.join(PATH_STATIC, 'test_example1.ort')) + # Verify the file indeed has no model + assert orso_data[0].info.data_source.sample.model is None + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + result = load_orso_model(orso_data) + + assert result is None + assert len(w) == 1 + assert 'does not contain a sample model definition' in str(w[0].message) diff --git a/tests/test_ort_file.py b/tests/test_ort_file.py index ff66c743..fe40e748 100644 --- a/tests/test_ort_file.py +++ b/tests/test_ort_file.py @@ -24,59 +24,57 @@ def make_pooch(base_url: str, registry: dict[str, str | None]) -> pooch.Pooch: """Make a Pooch object to download test data.""" return pooch.create( - path=pooch.os_cache("data"), - env="POOCH_DIR", + path=pooch.os_cache('data'), + env='POOCH_DIR', base_url=base_url, registry=registry, ) -@pytest.fixture(scope="module") +@pytest.fixture(scope='module') def data_registry(): return make_pooch( - base_url="https://pub-6c25ef91903d4301a3338bd53b370098.r2.dev", + base_url='https://pub-6c25ef91903d4301a3338bd53b370098.r2.dev', registry={ - "amor_reduced_iofq.ort": None, + 'amor_reduced_iofq.ort': None, }, ) -@pytest.fixture(scope="module") +@pytest.fixture(scope='module') def load_data(data_registry): - path = data_registry.fetch("amor_reduced_iofq.ort") - logging.info("Loading data from %s", path) + path = data_registry.fetch('amor_reduced_iofq.ort') + logging.info('Loading data from %s', path) data = load(path) return data -@pytest.fixture(scope="module") +@pytest.fixture(scope='module') def fit_model(load_data): data = load_data # Rescale data - reflectivity = data["data"]["R_0"].values + reflectivity = data['data']['R_0'].values scale_factor = 1 / np.max(reflectivity) - data["data"]["R_0"].values *= scale_factor + data['data']['R_0'].values *= scale_factor # Create a model for the sample - si = Material(sld=2.07, isld=0.0, name="Si") - sio2 = Material(sld=3.47, isld=0.0, name="SiO2") - d2o = Material(sld=6.33, isld=0.0, name="D2O") - dlipids = Material(sld=5.0, isld=0.0, name="DLipids") + si = Material(sld=2.07, isld=0.0, name='Si') + sio2 = Material(sld=3.47, isld=0.0, name='SiO2') + d2o = Material(sld=6.33, isld=0.0, name='D2O') + dlipids = Material(sld=5.0, isld=0.0, name='DLipids') - superphase = Layer(material=si, thickness=0, roughness=0, name="Si superphase") - sio2_layer = Layer(material=sio2, thickness=20, roughness=4, name="SiO2 layer") - dlipids_layer = Layer( - material=dlipids, thickness=40, roughness=4, name="DLipids layer" - ) - subphase = Layer(material=d2o, thickness=0, roughness=5, name="D2O subphase") + superphase = Layer(material=si, thickness=0, roughness=0, name='Si superphase') + sio2_layer = Layer(material=sio2, thickness=20, roughness=4, name='SiO2 layer') + dlipids_layer = Layer(material=dlipids, thickness=40, roughness=4, name='DLipids layer') + subphase = Layer(material=d2o, thickness=0, roughness=5, name='D2O subphase') multi_sample = Sample( Multilayer(superphase), Multilayer(sio2_layer), Multilayer(dlipids_layer), Multilayer(subphase), - name="Multilayer Structure", + name='Multilayer Structure', ) multi_layer_model = Model( @@ -84,7 +82,7 @@ def fit_model(load_data): scale=1, background=0.000001, resolution_function=PercentageFwhm(0), - name="Multilayer Model", + name='Multilayer Model', ) # Set the fitting parameters @@ -122,66 +120,60 @@ def fit_model(load_data): def test_read_reduced_data__check_structure(load_data): - data_keys = load_data["data"].keys() - coord_keys = load_data["coords"].keys() + data_keys = load_data['data'].keys() + coord_keys = load_data['coords'].keys() for key in data_keys: if key in coord_keys: - assert len(load_data["data"][key].values) == len( - load_data["coords"][key].values - ) + assert len(load_data['data'][key].values) == len(load_data['coords'][key].values) def test_validate_physical_data__r_values_non_negative(load_data): - for key in load_data["data"].keys(): - assert all(load_data["data"][key].values >= 0) + for key in load_data['data'].keys(): + assert all(load_data['data'][key].values >= 0) def test_validate_physical_data__r_values_finite(load_data): - for key in load_data["data"].keys(): - assert all(np.isfinite(load_data["data"][key].values)) + for key in load_data['data'].keys(): + assert all(np.isfinite(load_data['data'][key].values)) -@pytest.mark.skip("Currently no warning implemented") +@pytest.mark.skip('Currently no warning implemented') def test_validate_physical_data__r_values_ureal_positive(load_data): - a = load_data["data"]["R_0"].values - b = 1 + 2 * np.sqrt(load_data["data"]["R_0"].variances) + a = load_data['data']['R_0'].values + b = 1 + 2 * np.sqrt(load_data['data']['R_0'].variances) for val_a, val_b in zip(a, b): if val_a > val_b: pytest.warns( - UserWarning, - reason=f"Reflectivity value {val_a} is unphysically large compared to its uncertainty {val_b}" + UserWarning, reason=f'Reflectivity value {val_a} is unphysically large compared to its uncertainty {val_b}' ) - assert all( - load_data["data"]["R_0"].values - <= 1 + 2 * np.sqrt(load_data["data"]["R_0"].variances) - ) + assert all(load_data['data']['R_0'].values <= 1 + 2 * np.sqrt(load_data['data']['R_0'].variances)) def test_validate_physical_data__q_values_non_negative(load_data): - for key in load_data["coords"].keys(): - assert all(load_data["coords"][key].values >= 0) + for key in load_data['coords'].keys(): + assert all(load_data['coords'][key].values >= 0) def test_validate_physical_data__q_values_ureal_positive(load_data): - for key in load_data["coords"].keys(): + for key in load_data['coords'].keys(): # Reflectometry data is usually with the range of 0-5, # so 10 is a safe upper limit - assert all(load_data["coords"][key].values < 10) + assert all(load_data['coords'][key].values < 10) def test_validate_physical_data__q_values_finite(load_data): - for key in load_data["coords"].keys(): - assert all(np.isfinite(load_data["coords"][key].values < 10)) + for key in load_data['coords'].keys(): + assert all(np.isfinite(load_data['coords'][key].values < 10)) -@pytest.mark.skip("Currently no meta data to check") +@pytest.mark.skip('Currently no meta data to check') def test_validate_meta_data__required_meta_data() -> None: - pytest.fail(reason="Currently no meta data to check") + pytest.fail(reason='Currently no meta data to check') def test_analyze_reduced_data__fit_model_success(fit_model): - assert fit_model["success"] is True + assert fit_model['success'] is True def test_analyze_reduced_data__fit_model_reasonable(fit_model): - assert fit_model["reduced_chi"] < 0.01 + assert fit_model['reduced_chi'] < 0.01 diff --git a/tests/test_project.py b/tests/test_project.py index b86210ec..40febe90 100644 --- a/tests/test_project.py +++ b/tests/test_project.py @@ -4,6 +4,7 @@ from unittest.mock import MagicMock import numpy as np +import pytest from easyscience import global_object from easyscience.fitting import AvailableMinimizers from easyscience.variable import Parameter @@ -17,8 +18,11 @@ from easyreflectometry.model import PercentageFwhm from easyreflectometry.model import Pointwise from easyreflectometry.project import Project +from easyreflectometry.sample import Layer from easyreflectometry.sample import Material from easyreflectometry.sample import MaterialCollection +from easyreflectometry.sample import Multilayer +from easyreflectometry.sample import Sample PATH_STATIC = os.path.join(os.path.dirname(easyreflectometry.__file__), '..', '..', 'tests', '_static') @@ -576,6 +580,7 @@ def test_create(self, tmp_path): def test_load_experiment(self): # When + global_object.map._clear() project = Project() model_5 = Model() project.models = ModelCollection(Model(), Model(), Model(), Model(), Model(), model_5) @@ -587,13 +592,50 @@ def test_load_experiment(self): # Expect assert list(project.experiments.keys()) == [5] assert isinstance(project.experiments[5], DataSet1D) - assert project.experiments[5].name == 'Experiment 5' + assert project.experiments[5].name == 'Example data file from refnx docs' assert project.experiments[5].model == model_5 assert isinstance(project.models[5].resolution_function, Pointwise) assert isinstance(project.models[4].resolution_function, PercentageFwhm) + def test_load_experiment_sets_resolution_function_pointwise_when_xe_present(self, tmp_path): + # When + global_object.map._clear() + project = Project() + project.models = ModelCollection(Model()) + + # Create a simple 4-column data file (x, y, e, xe) + fpath = tmp_path / 'four_col.txt' + fpath.write_text('# test data\n0.01 1e-5 1e-6 1e-4\n0.02 2e-5 1e-6 1e-4\n') + + # Then + project.load_experiment_for_model_at_index(str(fpath)) + + # Expect Pointwise because xe (q-resolution) is present + from easyreflectometry.model.resolution_functions import Pointwise + + assert isinstance(project.models[0].resolution_function, Pointwise) + + def test_load_experiment_sets_linearspline_when_only_ye_present(self, tmp_path): + # When + global_object.map._clear() + project = Project() + project.models = ModelCollection(Model()) + + # Create a simple 3-column data file (x, y, e) + fpath = tmp_path / 'three_col.txt' + fpath.write_text('# test data\n0.01 1e-5 1e-6\n0.02 2e-5 1e-6\n') + + # Then + project.load_experiment_for_model_at_index(str(fpath)) + + # Expect LinearSpline because only ye (fwhm from ye) is available + from easyreflectometry.model.resolution_functions import LinearSpline + + assert isinstance(project.models[0].resolution_function, LinearSpline) + def test_experimental_data_at_index(self): # When + global_object.map._clear() project = Project() project.models = ModelCollection(Model()) fpath = os.path.join(PATH_STATIC, 'example.ort') @@ -603,7 +645,7 @@ def test_experimental_data_at_index(self): data = project.experimental_data_for_model_at_index() # Expect - assert data.name == 'Experiment 0' + assert data.name == 'Example data file from refnx docs' assert data.is_experiment assert isinstance(data, DataSet1D) assert len(data.x) == 408 @@ -613,6 +655,7 @@ def test_experimental_data_at_index(self): def test_q(self): # When + global_object.map._clear() project = Project() # Then @@ -688,15 +731,376 @@ def test_current_experiment_index_setter_out_of_range(self): project._experiments[0] = DataSet1D(name='exp0', x=[], y=[], ye=[], xe=[], model=None) # Negative index should raise - try: + with pytest.raises(ValueError): project.current_experiment_index = -1 - assert False, 'Expected ValueError for negative index' - except ValueError: - pass # Index >= len(_experiments) should raise - try: + with pytest.raises(ValueError): project.current_experiment_index = 1 - assert False, 'Expected ValueError for out-of-range index' - except ValueError: - pass + + def test_get_materials_from_model(self): + # When + global_object.map._clear() + project = Project() + material_1 = Material(sld=2.07, isld=0.0, name='Material 1') + material_2 = Material(sld=3.47, isld=0.0, name='Material 2') + material_3 = Material(sld=6.36, isld=0.0, name='Material 3') + + layer_1 = Layer(material=material_1, thickness=10, roughness=0, name='Layer 1') + layer_2 = Layer(material=material_2, thickness=20, roughness=1, name='Layer 2') + layer_3 = Layer(material=material_3, thickness=0, roughness=2, name='Layer 3') + + sample = Sample(Multilayer([layer_1, layer_2]), Multilayer([layer_3])) + model = Model(sample=sample) + + # Then + materials = project._get_materials_from_model(model) + + # Expect + assert len(materials) == 3 + assert materials[0] == material_1 + assert materials[1] == material_2 + assert materials[2] == material_3 + + def test_get_materials_from_model_duplicate_materials(self): + # When + global_object.map._clear() + project = Project() + # Use the same material in multiple layers + shared_material = Material(sld=2.07, isld=0.0, name='Shared Material') + material_2 = Material(sld=3.47, isld=0.0, name='Material 2') + + layer_1 = Layer(material=shared_material, thickness=10, roughness=0, name='Layer 1') + layer_2 = Layer(material=material_2, thickness=20, roughness=1, name='Layer 2') + layer_3 = Layer(material=shared_material, thickness=30, roughness=2, name='Layer 3') + + sample = Sample(Multilayer([layer_1, layer_2, layer_3])) + model = Model(sample=sample) + + # Then + materials = project._get_materials_from_model(model) + + # Expect - should only include unique materials + assert len(materials) == 2 + assert materials[0] == shared_material + assert materials[1] == material_2 + + def test_add_sample_from_orso(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + + initial_model_count = len(project._models) + initial_material_count = len(project._materials) + + material_1 = Material(sld=4.0, isld=0.0, name='New Material 1') + material_2 = Material(sld=5.0, isld=0.0, name='New Material 2') + layer_1 = Layer(material=material_1, thickness=50, roughness=1, name='New Layer 1') + layer_2 = Layer(material=material_2, thickness=100, roughness=2, name='New Layer 2') + new_sample = Sample(Multilayer([layer_1, layer_2])) + + # Then + project.add_sample_from_orso(new_sample) + + # Expect + assert len(project._models) == initial_model_count + 1 + assert project._models[-1].sample == new_sample + # The interface should be set by add_sample_from_orso + assert project._models[-1].interface == project._calculator + assert len(project._materials) == initial_material_count + 2 + assert material_1 in project._materials + assert material_2 in project._materials + assert project.current_model_index == len(project._models) - 1 + + def test_add_sample_from_orso_multiple_additions(self): + # When + global_object.map._clear() + project = Project() + + material_1 = Material(sld=2.0, isld=0.0, name='Material A') + layer_1 = Layer(material=material_1, thickness=10, roughness=0, name='Layer A') + sample_1 = Sample(Multilayer([layer_1])) + + material_2 = Material(sld=3.0, isld=0.0, name='Material B') + layer_2 = Layer(material=material_2, thickness=20, roughness=1, name='Layer B') + sample_2 = Sample(Multilayer([layer_2])) + + # Then + project.add_sample_from_orso(sample_1) + project.add_sample_from_orso(sample_2) + + # Expect + assert len(project._models) == 2 + assert project._models[0].sample == sample_1 + assert project._models[1].sample == sample_2 + assert len(project._materials) == 2 + assert material_1 in project._materials + assert material_2 in project._materials + assert project.current_model_index == 1 + + def test_add_sample_from_orso_with_shared_materials(self): + # When + global_object.map._clear() + project = Project() + + # Create first sample with a material + shared_material = Material(sld=2.0, isld=0.0, name='Shared Material') + layer_1 = Layer(material=shared_material, thickness=10, roughness=0, name='Layer 1') + sample_1 = Sample(Multilayer([layer_1])) + project.add_sample_from_orso(sample_1) + + initial_material_count = len(project._materials) + + # Create second sample using the same material + layer_2 = Layer(material=shared_material, thickness=20, roughness=1, name='Layer 2') + sample_2 = Sample(Multilayer([layer_2])) + + # Then + project.add_sample_from_orso(sample_2) + + # Expect - shared material should not be duplicated + assert len(project._models) == 2 + # The shared material instance is already in the collection, so count should stay the same + assert len(project._materials) == initial_material_count + + def test_replace_models_from_orso(self): + """Test that replace_models_from_orso replaces all existing models with a single new model.""" + # When + global_object.map._clear() + project = Project() + project.default_model() + + # Add some models to start with + material_1 = Material(sld=2.0, isld=0.0, name='Material 1') + layer_1 = Layer(material=material_1, thickness=10, roughness=0, name='Layer 1') + sample_1 = Sample(Multilayer([layer_1])) + project.add_sample_from_orso(sample_1) + + material_2 = Material(sld=3.0, isld=0.0, name='Material 2') + layer_2 = Layer(material=material_2, thickness=20, roughness=1, name='Layer 2') + sample_2 = Sample(Multilayer([layer_2])) + project.add_sample_from_orso(sample_2) + + # Verify we have multiple models + assert len(project._models) > 1 + len(project._models) + + # Create a new sample to replace all existing models + new_material = Material(sld=5.0, isld=0.5, name='New Material') + new_layer = Layer(material=new_material, thickness=50, roughness=2, name='New Layer') + new_sample = Sample(Multilayer([new_layer])) + + # Then - replace all models with the new sample + project.replace_models_from_orso(new_sample) + + # Expect - only one model should remain + assert len(project._models) == 1 + assert project._models[0].sample == new_sample + # The interface should be set + assert project._models[0].interface == project._calculator + # Only the new material should be in the materials collection + assert len(project._materials) == 1 + assert new_material in project._materials + # Old materials should not be in the collection + assert material_1 not in project._materials + assert material_2 not in project._materials + # Current model index should be reset to 0 + assert project.current_model_index == 0 + + def test_is_default_model_true(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + + # Then Expect + assert project.is_default_model(0) is True + + def test_is_default_model_false_non_default_sample(self): + # When + global_object.map._clear() + project = Project() + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project.models = ModelCollection(model) + + # Then Expect + assert project.is_default_model(0) is False + + def test_is_default_model_index_out_of_range(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + + # Then Expect + assert project.is_default_model(-1) is False + assert project.is_default_model(1) is False + assert project.is_default_model(100) is False + + def test_is_default_model_multiple_models(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a custom model + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + + # Then Expect + assert project.is_default_model(0) is True + assert project.is_default_model(1) is False + + def test_remove_model_at_index(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a second model + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + assert len(project._models) == 2 + + # Then + project.remove_model_at_index(0) + + # Expect + assert len(project._models) == 1 + assert project._models[0].sample[0].name == 'Custom Assembly' + + def test_remove_model_at_index_adjusts_current_index(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a second model + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + project._current_model_index = 1 + project._current_assembly_index = 1 + project._current_layer_index = 1 + + # Then + project.remove_model_at_index(0) + + # Expect - current_model_index should be adjusted + assert project._current_model_index == 0 + assert project._current_assembly_index == 0 + assert project._current_layer_index == 0 + + def test_remove_model_at_index_resets_indices_when_at_end(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a second model + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + project._current_model_index = 1 + + # Then - remove the model at the current index + project.remove_model_at_index(1) + + # Expect - current_model_index should be clamped to valid range + assert project._current_model_index == 0 + assert project._current_assembly_index == 0 + assert project._current_layer_index == 0 + + def test_remove_model_at_index_removes_experiment_at_same_index(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a second model + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + # Add experiment linked to model 0 + experiment = DataSet1D( + name='exp0', x=[0.01, 0.02], y=[1.0, 0.5], ye=[0.1, 0.1], xe=[0.001, 0.001], model=project._models[0] + ) + project._experiments[0] = experiment + + # Then + project.remove_model_at_index(0) + + # Expect - experiment mapped to the removed model index is removed + assert 0 not in project._experiments + + def test_remove_model_at_index_reindexes_experiments_above_removed_index(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + + # Add two more models (total = 3) + material_1 = Material(sld=4.0, isld=0.0, name='Custom Material 1') + layer_1 = Layer(material=material_1, thickness=50, roughness=1, name='Custom Layer 1') + model_1 = Model(sample=Sample(Multilayer([layer_1], name='Custom Assembly 1'))) + project._models.append(model_1) + + material_2 = Material(sld=5.0, isld=0.0, name='Custom Material 2') + layer_2 = Layer(material=material_2, thickness=60, roughness=2, name='Custom Layer 2') + model_2 = Model(sample=Sample(Multilayer([layer_2], name='Custom Assembly 2'))) + project._models.append(model_2) + + # Add experiments for all model indices 0, 1, 2 + project._experiments[0] = DataSet1D(name='exp0', x=[0.01], y=[1.0], ye=[0.1], xe=[0.001], model=project._models[0]) + project._experiments[1] = DataSet1D(name='exp1', x=[0.02], y=[0.9], ye=[0.1], xe=[0.001], model=project._models[1]) + project._experiments[2] = DataSet1D(name='exp2', x=[0.03], y=[0.8], ye=[0.1], xe=[0.001], model=project._models[2]) + + # Then - remove middle model + project.remove_model_at_index(1) + + # Expect - middle experiment removed and upper one shifted down + assert set(project._experiments.keys()) == {0, 1} + assert project._experiments[0].name == 'exp0' + assert project._experiments[1].name == 'exp2' + + def test_remove_model_at_index_raises_for_last_model(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + assert len(project._models) == 1 + + # Then Expect + with pytest.raises(ValueError, match='Cannot remove the last model'): + project.remove_model_at_index(0) + + def test_remove_model_at_index_raises_for_invalid_index(self): + # When + global_object.map._clear() + project = Project() + project.default_model() + # Add a second model so we have 2 + material = Material(sld=4.0, isld=0.0, name='Custom Material') + layer = Layer(material=material, thickness=50, roughness=1, name='Custom Layer') + sample = Sample(Multilayer([layer], name='Custom Assembly')) + model = Model(sample=sample) + project._models.append(model) + + # Then Expect - negative index + with pytest.raises(IndexError, match='out of range'): + project.remove_model_at_index(-1) + + # Then Expect - index >= len + with pytest.raises(IndexError, match='out of range'): + project.remove_model_at_index(2) From b65486cb8aad2a72142cf46e30958ad452e09ec2 Mon Sep 17 00:00:00 2001 From: Piotr Rozyczko Date: Fri, 27 Feb 2026 20:40:18 +0100 Subject: [PATCH 8/9] Improved sample handling (#293) * fix for appending samples from loaded models * added unit tests * use most recent version of refl1d * update docs to include add_sample_from_orso * temporarily revert to fwhm * minor fixes after CR + ruff format * back to pointwise * proper conversion from stack to layers. * added method * temporarily switch off pointwise * fixed SLD read. * added handling for orso names (model and experiment) * minor fix for names * CR review fixes * model reindexing issue fixed. Project rst added * path for SLD-less orso files * fix variances sent to bumps. Make FWHM default (temporarily) * fixed test to correspond to new resolution * added logs for changing minimizers * reload constraints where necessary * expose chi2 properly * PR review #1 * Fix the zero variances issue * allow 3.13 * added some tests --- .copier-answers.yml | 2 +- .github/workflows/python-ci.yml | 2 +- .github/workflows/python-package.yml | 2 +- CONTRIBUTING.rst | 2 +- pyproject.toml | 6 +- src/easyreflectometry/fitting.py | 50 ++++++++- src/easyreflectometry/project.py | 39 ++++--- tests/summary/test_summary.py | 2 +- tests/test_fitting.py | 152 +++++++++++++++++++++++++++ tests/test_ort_file.py | 33 +++--- tests/test_project.py | 16 +-- 11 files changed, 258 insertions(+), 48 deletions(-) diff --git a/.copier-answers.yml b/.copier-answers.yml index 0653ab1c..867a00e8 100644 --- a/.copier-answers.yml +++ b/.copier-answers.yml @@ -2,7 +2,7 @@ _commit: 8bdcedc _src_path: gh:/EasyScience/EasyProjectTemplate description: A reflectometry python package built on the EasyScience framework. -max_python: '3.12' +max_python: '3.13' min_python: '3.9' orgname: EasyScience packagename: easyreflectometry diff --git a/.github/workflows/python-ci.yml b/.github/workflows/python-ci.yml index cea64fc2..90967cd9 100644 --- a/.github/workflows/python-ci.yml +++ b/.github/workflows/python-ci.yml @@ -28,7 +28,7 @@ jobs: strategy: max-parallel: 4 matrix: - python-version: ['3.11', '3.12'] + python-version: ['3.11', '3.12', '3.13'] os: [ubuntu-latest, macos-latest, windows-2022] runs-on: ${{ matrix.os }} diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index c14850d9..3a4d1036 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -18,7 +18,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ['3.11','3.12'] + python-version: ['3.11','3.12','3.13'] if: "!contains(github.event.head_commit.message, '[ci skip]')" steps: diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index e4a56d06..078ad7a9 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -102,7 +102,7 @@ Before you submit a pull request, check that it meets these guidelines: 2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring, and add the feature to the list in README.md. -3. The pull request should work for Python, 3.11 and 3.12, and for PyPy. Check +3. The pull request should work for Python, 3.11, 3.12, and 3.13, and for PyPy. Check https://travis-ci.com/easyScience/EasyReflectometryLib/pull_requests and make sure that the tests pass for all supported Python versions. diff --git a/pyproject.toml b/pyproject.toml index 630422ca..2c442dfe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,10 +23,11 @@ classifiers = [ "Programming Language :: Python :: 3 :: Only", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", "Development Status :: 3 - Alpha" ] -requires-python = ">=3.11,<3.13" +requires-python = ">=3.11,<3.14" dependencies = [ "easyscience @ git+https://github.com/easyscience/corelib.git@develop", @@ -134,11 +135,12 @@ force-single-line = true legacy_tox_ini = """ [tox] isolated_build = True -envlist = py{3.11,3.12} +envlist = py{3.11,3.12,3.13} [gh-actions] python = 3.11: py311 3.12: py312 + 3.13: py313 [gh-actions:env] PLATFORM = ubuntu-latest: linux diff --git a/src/easyreflectometry/fitting.py b/src/easyreflectometry/fitting.py index a33e6570..86df41e3 100644 --- a/src/easyreflectometry/fitting.py +++ b/src/easyreflectometry/fitting.py @@ -31,6 +31,7 @@ def wrapped(*args, **kwargs): self._fit_func = [func_wrapper(m.interface.fit_func, m.unique_name) for m in args] self._models = args self.easy_science_multi_fitter = EasyScienceMultiFitter(args, self._fit_func) + self._fit_results: list[FitResults] | None = None def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup: """ @@ -75,6 +76,7 @@ def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup: dy.append(1 / np.sqrt(variances_masked)) result = self.easy_science_multi_fitter.fit(x, y, weights=dy) + self._fit_results = result new_data = data.copy() for i, _ in enumerate(result): id = refl_nums[i] @@ -99,7 +101,53 @@ def fit_single_data_set_1d(self, data: DataSet1D) -> FitResults: :param data: DataGroup to be fitted to and populated :param method: Optimisation method """ - return self.easy_science_multi_fitter.fit(x=[data.x], y=[data.y], weights=[data.ye])[0] + x_vals = np.asarray(data.x) + y_vals = np.asarray(data.y) + variances = np.asarray(data.ye) + + zero_variance_mask = variances == 0.0 + num_zero_variance = int(np.sum(zero_variance_mask)) + + if num_zero_variance > 0: + warnings.warn( + f'Masked {num_zero_variance} data point(s) in single-dataset fit due to zero variance during fitting.', + UserWarning, + ) + + valid_mask = ~zero_variance_mask + if not np.any(valid_mask): + raise ValueError('Cannot fit single dataset: all points have zero variance.') + + x_vals_masked = x_vals[valid_mask] + y_vals_masked = y_vals[valid_mask] + variances_masked = variances[valid_mask] + + weights = 1.0 / np.sqrt(variances_masked) + result = self.easy_science_multi_fitter.fit(x=[x_vals_masked], y=[y_vals_masked], weights=[weights])[0] + self._fit_results = [result] + return result + + @property + def chi2(self) -> float | None: + """Total chi-squared across all fitted datasets, or None if no fit has been performed.""" + if self._fit_results is None: + return None + return sum(r.chi2 for r in self._fit_results) + + @property + def reduced_chi(self) -> float | None: + """Reduced chi-squared from the most recent fit, or None if no fit has been performed.""" + if self._fit_results is None: + return None + total_chi2 = sum(r.chi2 for r in self._fit_results) + total_points = sum(np.size(r.x) for r in self._fit_results) + n_params = self._fit_results[0].n_pars + total_dof = total_points - n_params + + if total_dof <= 0: + return None + + return total_chi2 / total_dof def switch_minimizer(self, minimizer: AvailableMinimizers) -> None: """ diff --git a/src/easyreflectometry/project.py b/src/easyreflectometry/project.py index 13b4fdda..8c7b6a73 100644 --- a/src/easyreflectometry/project.py +++ b/src/easyreflectometry/project.py @@ -1,5 +1,6 @@ import datetime import json +import logging import os from pathlib import Path from typing import Dict @@ -10,8 +11,8 @@ import numpy as np from easyscience import global_object from easyscience.fitting import AvailableMinimizers -from easyscience.fitting.fitter import DEFAULT_MINIMIZER from easyscience.variable import Parameter +from easyscience.variable.parameter_dependency_resolver import resolve_all_parameter_dependencies from scipp import DataGroup from easyreflectometry.calculators import CalculatorFactory @@ -20,11 +21,9 @@ from easyreflectometry.data.measurement import extract_orso_title from easyreflectometry.data.measurement import load_data_from_orso_file from easyreflectometry.fitting import MultiFitter -from easyreflectometry.model import LinearSpline from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm -from easyreflectometry.model import Pointwise from easyreflectometry.sample import Layer from easyreflectometry.sample import Material from easyreflectometry.sample import MaterialCollection @@ -32,11 +31,13 @@ from easyreflectometry.sample import Sample from easyreflectometry.sample.collections.base_collection import BaseCollection +logger = logging.getLogger(__name__) + Q_MIN = 0.001 Q_MAX = 0.3 Q_RESOLUTION = 500 -DEFAULT_MINIZER = AvailableMinimizers.LMFit_leastsq +DEFAULT_MINIMIZER = AvailableMinimizers.LMFit_leastsq class Project: @@ -48,6 +49,7 @@ def __init__(self): self._calculator = CalculatorFactory() self._experiments: Dict[DataGroup] = {} self._fitter: MultiFitter = None + self._minimizer_selection: AvailableMinimizers = DEFAULT_MINIMIZER self._colors: list[str] = None self._report = None self._q_min: float = None @@ -207,9 +209,8 @@ def models(self, models: ModelCollection) -> None: def fitter(self) -> MultiFitter: if len(self._models): if (self._fitter is None) or (self._fitter_model_index != self._current_model_index): - minimizer = self.minimizer self._fitter = MultiFitter(self._models[self._current_model_index]) - self.minimizer = minimizer + self._fitter.easy_science_multi_fitter.switch_minimizer(self._minimizer_selection) self._fitter_model_index = self._current_model_index return self._fitter @@ -225,10 +226,14 @@ def calculator(self, calculator: str) -> None: def minimizer(self) -> AvailableMinimizers: if self._fitter is not None: return self._fitter.easy_science_multi_fitter.minimizer.enum - return DEFAULT_MINIMIZER + return self._minimizer_selection @minimizer.setter def minimizer(self, minimizer: AvailableMinimizers) -> None: + old_name = getattr(self._minimizer_selection, 'name', str(self._minimizer_selection)) + new_name = getattr(minimizer, 'name', str(minimizer)) + logger.info('Minimizer changed from %s to %s (fitter active: %s)', old_name, new_name, self._fitter is not None) + self._minimizer_selection = minimizer if self._fitter is not None: self._fitter.easy_science_multi_fitter.switch_minimizer(minimizer) @@ -386,21 +391,10 @@ def _apply_resolution_function( ) -> None: """Set the resolution function on *model* based on variance data in *experiment*. - Prefers Pointwise when q-resolution (xe) data is present, otherwise falls - back to LinearSpline when reflectivity error (ye) data is present. - :param experiment: The experiment whose variance data drives the choice. :param model: The model whose resolution function is set. """ - if sum(experiment.xe) != 0: - resolution_function = Pointwise(q_data_points=[experiment.x, experiment.y, experiment.xe]) - model.resolution_function = resolution_function - elif sum(experiment.ye) != 0: - resolution_function = LinearSpline( - q_data_points=experiment.x, - fwhm_values=np.sqrt(experiment.ye), - ) - model.resolution_function = resolution_function + model.resolution_function = PercentageFwhm(5.0) def load_new_experiment(self, path: Union[Path, str]) -> None: new_experiment = load_as_dataset(str(path)) @@ -603,6 +597,8 @@ def as_dict(self, include_materials_not_in_model=False): self._as_dict_add_experiments(project_dict) if self.fitter is not None: project_dict['fitter_minimizer'] = self.fitter.easy_science_multi_fitter.minimizer.name + elif self._minimizer_selection is not None: + project_dict['fitter_minimizer'] = self._minimizer_selection.name if self._calculator is not None: project_dict['calculator'] = self._calculator.current_interface_name if self._colors is not None: @@ -641,7 +637,7 @@ def from_dict(self, project_dict: dict): if 'materials_not_in_model' in keys: self._materials.extend(MaterialCollection.from_dict(project_dict['materials_not_in_model'])) if 'fitter_minimizer' in keys: - self.fitter.easy_science_multi_fitter.switch_minimizer(AvailableMinimizers[project_dict['fitter_minimizer']]) + self.minimizer = AvailableMinimizers[project_dict['fitter_minimizer']] else: self._fitter = None if 'experiments' in keys: @@ -649,6 +645,9 @@ def from_dict(self, project_dict: dict): else: self._experiments = {} + # Resolve any pending parameter dependencies (constraints) after all objects are loaded + resolve_all_parameter_dependencies(self) + def _from_dict_extract_experiments(self, project_dict: dict) -> Dict[int, DataSet1D]: experiments = {} for key in project_dict['experiments'].keys(): diff --git a/tests/summary/test_summary.py b/tests/summary/test_summary.py index 65de9ba3..3bbb186e 100644 --- a/tests/summary/test_summary.py +++ b/tests/summary/test_summary.py @@ -133,7 +133,7 @@ def test_experiments_section(self, project: Project) -> None: assert 'No. of data points' in html assert '408' in html assert 'Resolution function' in html - assert 'Pointwise' in html + assert 'PercentageFwhm' in html def test_experiments_section_percentage_fhwm(self, project: Project) -> None: # When diff --git a/tests/test_fitting.py b/tests/test_fitting.py index fdeba385..446f10b9 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -1,12 +1,15 @@ __author__ = 'github.com/arm61' import os +from unittest.mock import MagicMock +import numpy as np import pytest from easyscience.fitting.minimizers.factory import AvailableMinimizers import easyreflectometry from easyreflectometry.calculators import CalculatorFactory +from easyreflectometry.data import DataSet1D from easyreflectometry.data.measurement import load from easyreflectometry.fitting import MultiFitter from easyreflectometry.model import Model @@ -224,3 +227,152 @@ def test_fitting_with_manual_zero_variance(): assert 'R_0_model' in analysed.keys() assert 'SLD_0' in analysed.keys() assert 'success' in analysed.keys() + + +def test_fit_single_data_set_1d_masks_zero_variance_points(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + captured = {} + mock_result = MagicMock() + mock_result.chi2 = 1.0 + mock_result.n_pars = 1 + + def _fake_fit(*, x, y, weights): + captured['x'] = x + captured['y'] = y + captured['weights'] = weights + return [mock_result] + + fitter.easy_science_multi_fitter = MagicMock() + fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit) + + data = DataSet1D( + name='single_dataset', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.01, 0.0, 0.04]), + ) + + with pytest.warns(UserWarning, match='Masked 1 data point\(s\) in single-dataset fit'): + result = fitter.fit_single_data_set_1d(data) + + assert result is mock_result + assert np.allclose(captured['x'][0], np.array([0.01, 0.03])) + assert np.allclose(captured['y'][0], np.array([1.0, 0.6])) + assert np.allclose(captured['weights'][0], np.array([10.0, 5.0])) + + +def test_reduced_chi_uses_global_dof_across_fit_results(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + fit_result_1 = MagicMock() + fit_result_1.chi2 = 10.0 + fit_result_1.x = np.arange(5) + fit_result_1.n_pars = 3 + + fit_result_2 = MagicMock() + fit_result_2.chi2 = 14.0 + fit_result_2.x = np.arange(7) + fit_result_2.n_pars = 3 + + fitter._fit_results = [fit_result_1, fit_result_2] + + expected = (10.0 + 14.0) / ((5 + 7) - 3) + assert fitter.reduced_chi == pytest.approx(expected) + + +def test_fit_single_data_set_1d_all_zero_variance_raises(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + data = DataSet1D( + name='all_zero', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.0, 0.0, 0.0]), + ) + + with pytest.raises(ValueError, match='all points have zero variance'): + fitter.fit_single_data_set_1d(data) + + +def test_chi2_returns_none_before_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + assert fitter.chi2 is None + + +def test_chi2_returns_total_after_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + r1 = MagicMock() + r1.chi2 = 5.0 + r2 = MagicMock() + r2.chi2 = 3.0 + + fitter._fit_results = [r1, r2] + assert fitter.chi2 == pytest.approx(8.0) + + +def test_reduced_chi_returns_none_before_fit(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + assert fitter.reduced_chi is None + + +def test_reduced_chi_returns_none_when_dof_zero(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + r1 = MagicMock() + r1.chi2 = 5.0 + r1.x = np.arange(3) + r1.n_pars = 3 # total_points == n_params => dof == 0 + + fitter._fit_results = [r1] + assert fitter.reduced_chi is None + + +def test_fit_single_data_set_1d_no_zero_variance(): + model = Model() + model.interface = CalculatorFactory() + fitter = MultiFitter(model) + + captured = {} + mock_result = MagicMock() + mock_result.chi2 = 2.0 + mock_result.n_pars = 1 + + def _fake_fit(*, x, y, weights): + captured['x'] = x + captured['y'] = y + captured['weights'] = weights + return [mock_result] + + fitter.easy_science_multi_fitter = MagicMock() + fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit) + + data = DataSet1D( + name='no_zero', + x=np.array([0.01, 0.02, 0.03]), + y=np.array([1.0, 0.8, 0.6]), + ye=np.array([0.01, 0.04, 0.09]), + ) + + result = fitter.fit_single_data_set_1d(data) + + assert result is mock_result + assert np.allclose(captured['x'][0], np.array([0.01, 0.02, 0.03])) + assert np.allclose(captured['y'][0], np.array([1.0, 0.8, 0.6])) diff --git a/tests/test_ort_file.py b/tests/test_ort_file.py index fe40e748..c547b1f5 100644 --- a/tests/test_ort_file.py +++ b/tests/test_ort_file.py @@ -56,6 +56,7 @@ def fit_model(load_data): reflectivity = data['data']['R_0'].values scale_factor = 1 / np.max(reflectivity) data['data']['R_0'].values *= scale_factor + data['data']['R_0'].variances *= scale_factor**2 # Create a model for the sample @@ -81,22 +82,30 @@ def fit_model(load_data): sample=multi_sample, scale=1, background=0.000001, - resolution_function=PercentageFwhm(0), + resolution_function=PercentageFwhm(5), name='Multilayer Model', ) # Set the fitting parameters - sio2_layer.roughness.bounds = (3, 12) - sio2_layer.material.sld.bounds = (3.47, 5) - sio2_layer.thickness.bounds = (10, 30) - - subphase.material.sld.bounds = (6, 6.35) - dlipids_layer.thickness.bounds = (30, 60) - dlipids_layer.roughness.bounds = (3, 10) - dlipids_layer.material.sld.bounds = (4, 6) - multi_layer_model.scale.bounds = (0.8, 1.2) - multi_layer_model.background.bounds = (1e-6, 1e-3) + sio2_layer.roughness.min = 3 + sio2_layer.roughness.max = 12 + sio2_layer.material.sld.min = 3.47 + sio2_layer.material.sld.max = 5 + sio2_layer.thickness.min = 10 + sio2_layer.thickness.max = 30 + + subphase.material.sld.min = 6 + dlipids_layer.thickness.min = 30 + dlipids_layer.thickness.max = 60 + dlipids_layer.roughness.min = 3 + dlipids_layer.roughness.max = 10 + dlipids_layer.material.sld.min = 4 + dlipids_layer.material.sld.max = 6 + multi_layer_model.scale.min = 0.8 + multi_layer_model.scale.max = 1.2 + multi_layer_model.background.min = 1e-6 + multi_layer_model.background.max = 1e-3 sio2_layer.roughness.free = True sio2_layer.material.sld.free = True @@ -176,4 +185,4 @@ def test_analyze_reduced_data__fit_model_success(fit_model): def test_analyze_reduced_data__fit_model_reasonable(fit_model): - assert fit_model['reduced_chi'] < 0.01 + assert fit_model['reduced_chi'] < 6.0 diff --git a/tests/test_project.py b/tests/test_project.py index 40febe90..b93df068 100644 --- a/tests/test_project.py +++ b/tests/test_project.py @@ -16,7 +16,6 @@ from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm -from easyreflectometry.model import Pointwise from easyreflectometry.project import Project from easyreflectometry.sample import Layer from easyreflectometry.sample import Material @@ -350,6 +349,7 @@ def test_as_dict(self): keys.sort() assert keys == [ 'calculator', + 'fitter_minimizer', 'info', 'models', 'with_experiments', @@ -594,7 +594,7 @@ def test_load_experiment(self): assert isinstance(project.experiments[5], DataSet1D) assert project.experiments[5].name == 'Example data file from refnx docs' assert project.experiments[5].model == model_5 - assert isinstance(project.models[5].resolution_function, Pointwise) + assert isinstance(project.models[5].resolution_function, PercentageFwhm) assert isinstance(project.models[4].resolution_function, PercentageFwhm) def test_load_experiment_sets_resolution_function_pointwise_when_xe_present(self, tmp_path): @@ -610,10 +610,10 @@ def test_load_experiment_sets_resolution_function_pointwise_when_xe_present(self # Then project.load_experiment_for_model_at_index(str(fpath)) - # Expect Pointwise because xe (q-resolution) is present - from easyreflectometry.model.resolution_functions import Pointwise + # Resolution is always set to PercentageFwhm + from easyreflectometry.model.resolution_functions import PercentageFwhm - assert isinstance(project.models[0].resolution_function, Pointwise) + assert isinstance(project.models[0].resolution_function, PercentageFwhm) def test_load_experiment_sets_linearspline_when_only_ye_present(self, tmp_path): # When @@ -628,10 +628,10 @@ def test_load_experiment_sets_linearspline_when_only_ye_present(self, tmp_path): # Then project.load_experiment_for_model_at_index(str(fpath)) - # Expect LinearSpline because only ye (fwhm from ye) is available - from easyreflectometry.model.resolution_functions import LinearSpline + # Resolution is always set to PercentageFwhm + from easyreflectometry.model.resolution_functions import PercentageFwhm - assert isinstance(project.models[0].resolution_function, LinearSpline) + assert isinstance(project.models[0].resolution_function, PercentageFwhm) def test_experimental_data_at_index(self): # When From e6c06f8dd33653660c1f82be036162b895e7e83b Mon Sep 17 00:00:00 2001 From: Piotr Rozyczko Date: Fri, 27 Feb 2026 20:50:09 +0100 Subject: [PATCH 9/9] ruff --- src/easyreflectometry/project.py | 1 - tests/sample/assemblies/test_bilayer.py | 2 ++ 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/easyreflectometry/project.py b/src/easyreflectometry/project.py index a9a12fa4..8c7b6a73 100644 --- a/src/easyreflectometry/project.py +++ b/src/easyreflectometry/project.py @@ -21,7 +21,6 @@ from easyreflectometry.data.measurement import extract_orso_title from easyreflectometry.data.measurement import load_data_from_orso_file from easyreflectometry.fitting import MultiFitter -from easyreflectometry.model import LinearSpline from easyreflectometry.model import Model from easyreflectometry.model import ModelCollection from easyreflectometry.model import PercentageFwhm diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py index 1c577c94..b96e6a78 100644 --- a/tests/sample/assemblies/test_bilayer.py +++ b/tests/sample/assemblies/test_bilayer.py @@ -17,12 +17,14 @@ class TestBilayer: def setup_method(self): from easyscience import global_object + # Clear the global object map to prevent name collisions # Accessing private _clear method as Map doesn't expose a public clear if hasattr(global_object.map, 'clear'): global_object.map.clear() elif hasattr(global_object.map, '_clear'): global_object.map._clear() + def test_default(self): """Test default bilayer creation with expected structure.""" p = Bilayer()