diff --git a/docs/docs/tutorials/basic/material_library.md b/docs/docs/tutorials/basic/material_library.md index 78bdfcf1..655263a6 100644 --- a/docs/docs/tutorials/basic/material_library.md +++ b/docs/docs/tutorials/basic/material_library.md @@ -49,6 +49,34 @@ si = MaterialDensity(chemical_structure=chemical_structure, density=2.65, name=' The density should be in units of grams per cubic centimeter and the scattering length is calculated from `'SiO2'`. +By default the `sld` and `isld` of a `MaterialDensity` are _dependent_ +parameters, recomputed from the density, the formula's scattering length +and its molecular weight whenever any of those change - so `density` is +the parameter to vary in a fit, and assigning to `sld` directly is not +possible. This coupling can be switched off per material with the +`sld_coupled` property: + +```python +si.sld_coupled = False # sld/isld become independent, keep their values +si.sld.fixed = False # ...and can now be fitted directly +``` + +While decoupled, changes to `density` (or the formula) no longer +propagate to the SLD, and the density, molecular weight and scattering +length no longer affect the reflectivity. Setting `sld_coupled = True` +restores the dependency and **recalculates** `sld`/`isld` from the +current formula and density, discarding any manually set or fitted +values. The coupling state - and, when decoupled, the manual SLD +values - survive serialization (`as_dict`/`from_dict`); dictionaries +from before this feature restore as coupled. + +Assigning a new `chemical_structure` updates both the scattering length +and the molecular weight; a formula that does not parse to at least one +known atom raises `ValueError` and leaves the material unchanged. + +Note that `molecular_weight` is a read-only descriptor, not a fit +parameter: it is fully determined by the formula. + ## MaterialSolvated Sometimes it is desirable to have a layer that consists of a material diff --git a/src/easyreflectometry/sample/elements/materials/material_density.py b/src/easyreflectometry/sample/elements/materials/material_density.py index 5be1cf9d..1d6d046e 100644 --- a/src/easyreflectometry/sample/elements/materials/material_density.py +++ b/src/easyreflectometry/sample/elements/materials/material_density.py @@ -6,6 +6,7 @@ import numpy as np from easyscience import global_object +from easyscience.variable import DescriptorNumber from easyscience.variable import Parameter from easyreflectometry.special.calculations import density_to_sld @@ -27,20 +28,53 @@ 'max': np.inf, 'fixed': True, }, + # A DescriptorNumber, not a Parameter: the molecular weight is a constant + # of the chemical formula (recomputed whenever the formula changes) and + # must never enter a fit — it is fully degenerate with density in the + # derived SLD (only the ratio density/molecular_weight is observable). 'molecular_weight': { 'description': 'The molecular weight of a material.', 'url': 'https://en.wikipedia.org/wiki/Molecular_mass', 'value': 28.02, 'unit': 'g / mole', - 'min': -np.inf, - 'max': np.inf, - 'fixed': True, }, } DEFAULTS.update(MATERIAL_DEFAULTS) class MaterialDensity(Material): + """A material defined by chemical formula and mass density. + + The scattering length density is derived rather than set: from the + formula, the coherent neutron scattering length ``b`` (real and + imaginary parts, tabulated per isotope) and the molecular weight ``M`` + are computed, and ``sld``/``isld`` are wired as *dependent* parameters + + sld = N_A * density * b / M + + so ``density`` is the natural fit parameter and edits to the density or + the formula propagate to the SLD automatically. + + The coupling can be switched off per material via :attr:`sld_coupled`: + when ``False``, ``sld``/``isld`` are independent parameters that can be + set and fitted directly, while ``density``, ``molecular_weight`` and the + scattering lengths stop affecting anything until the coupling is + restored. Restoring it (``sld_coupled = True``) recomputes the SLDs from + the current formula and density, discarding manually set values. The + state round-trips through ``as_dict``/``from_dict``, including the + manual SLD values of a decoupled material; dictionaries from before + this feature deserialize as coupled. + + Assigning :attr:`chemical_structure` updates the scattering lengths and + the molecular weight together; an invalid formula raises ``ValueError`` + and leaves the material unchanged. + + :attr:`molecular_weight` is a read-only ``DescriptorNumber``, never a fit + parameter: it is fully determined by the formula, and freeing it alongside + density would make the fit degenerate (only ``density / molecular_weight`` + enters the derived SLD). + """ + def __init__( self, chemical_structure: Union[str, None] = None, @@ -79,11 +113,13 @@ def __init__( scattering_length = neutron_scattering_length(chemical_structure) - mw = get_as_parameter( + mw = DescriptorNumber( name='molecular_weight', value=molecular_weight(chemical_structure), - default_dict=DEFAULTS, - unique_name_prefix=f'{unique_name}_Mw', + unit=DEFAULTS['molecular_weight']['unit'], + description=DEFAULTS['molecular_weight']['description'], + url=DEFAULTS['molecular_weight']['url'], + unique_name=global_object.generate_unique_name(f'{unique_name}_Mw'), ) scattering_length_real = get_as_parameter( name='scattering_length_real', @@ -157,6 +193,41 @@ def _setup_sld_constraints(self) -> None: }, ) + @property + def sld_coupled(self) -> bool: + """Whether ``sld``/``isld`` are derived from formula & density (True, + the default) or independent, directly editable/fittable parameters + (False). The dependency state itself is the source of truth.""" + return not self._sld.independent + + @sld_coupled.setter + def sld_coupled(self, couple: bool) -> None: + if couple == self.sld_coupled: + return + if couple: + # Recomputes sld/isld from the current density/scattering + # length/molecular weight — manually set values are discarded. + self._setup_sld_constraints() + else: + # make_independent raises on an already-independent parameter, + # so guard each individually. Values are kept. + for parameter in (self._sld, self._isld): + if not parameter.independent: + parameter.make_independent() + + def _convert_to_dict(self, d: dict, serializer, skip: Optional[list] = None, **kwargs) -> dict: + """Serializer hook (see ``SerializerBase._convert_to_dict``). + + ``sld``/``isld`` are not constructor arguments, so the argspec-driven + encoder never persists them; in the decoupled state their manually + entered or fitted values would be lost on save/load without this. + """ + d['sld_coupled'] = self.sld_coupled + if not self.sld_coupled: + d['sld'] = self._sld.value + d['isld'] = self._isld.value + return d + @classmethod def from_dict(cls, obj_dict: dict) -> 'MaterialDensity': """Re-attach sld/isld dependencies after deserialization. @@ -166,9 +237,28 @@ def from_dict(cls, obj_dict: dict) -> 'MaterialDensity': the constraint graph built in `__init__` still references the temporary Parameter created from the float kwarg. Rebuild here so `q.density = X` propagates to the derived SLDs. + + The keys written by ``_convert_to_dict`` are not constructor + arguments and must be removed before the parent's ``cls(**data)`` + call; they are then used to restore the coupling state. A dict + without them (pre-feature project files) restores as coupled. """ + obj_dict = dict(obj_dict) + sld_coupled = obj_dict.pop('sld_coupled', True) + manual_sld = obj_dict.pop('sld', None) + manual_isld = obj_dict.pop('isld', None) + instance = super().from_dict(obj_dict) - instance._setup_sld_constraints() + if sld_coupled: + instance._setup_sld_constraints() + else: + # __init__ wired the dependencies; undo them and restore the + # saved manual values. + instance.sld_coupled = False + if manual_sld is not None: + instance._sld.value = manual_sld + if manual_isld is not None: + instance._isld.value = manual_isld return instance @property @@ -185,10 +275,21 @@ def chemical_structure(self, structure_string: str) -> None: structure_string : str String that defines the chemical structure. """ - self._chemical_structure = structure_string + # Derive everything before mutating any state: an invalid formula + # must leave the material fully unchanged. periodictable parses + # garbage to an *empty* formula (b=0, mw=0) instead of raising, and + # mw=0 would put a division by zero into the sld dependency. scattering_length = neutron_scattering_length(structure_string) + # The molecular weight enters the sld dependency alongside the + # scattering length; leaving it at the old formula's value would make + # the derived sld a mix of two formulas. + mw = molecular_weight(structure_string) + if not mw: + raise ValueError(f'Invalid chemical formula: {structure_string!r}') + self._chemical_structure = structure_string self._scattering_length_real.value = scattering_length.real self._scattering_length_imag.value = scattering_length.imag + self._molecular_weight.value = mw @property def density(self) -> Parameter: @@ -199,7 +300,10 @@ def density(self, value: float) -> None: self._density.value = value @property - def molecular_weight(self) -> Parameter: + def molecular_weight(self) -> DescriptorNumber: + """The molecular weight of the formula. A read-only descriptor, not a + fittable parameter: it is a constant of the chemical formula and is + recomputed whenever :attr:`chemical_structure` is assigned.""" return self._molecular_weight @property diff --git a/tests/sample/elements/materials/test_material_density.py b/tests/sample/elements/materials/test_material_density.py index 0b424cad..11a5677e 100644 --- a/tests/sample/elements/materials/test_material_density.py +++ b/tests/sample/elements/materials/test_material_density.py @@ -43,8 +43,11 @@ def test_chemical_structure_change(self): assert p.chemical_structure == 'Co' p.chemical_structure = 'B' assert p.density.value == 8.9 - assert_almost_equal(p.sld.value, 4.82010833570636) - assert_almost_equal(p.isld.value, -0.19098540517806603) + # The setter updates the molecular weight along with the scattering + # lengths; the derived sld reflects boron's mw, not cobalt's. + assert_almost_equal(p.molecular_weight.value, 10.81) + assert_almost_equal(p.sld.value, 26.277925961998147) + assert_almost_equal(p.isld.value, -1.0412008400037) assert p.chemical_structure == 'B' def test_dict_repr(self): @@ -65,6 +68,99 @@ def test_dict_round_trip(self): assert sorted(p.as_dict()) == sorted(q.as_dict()) + def test_chemical_structure_invalid_formula_leaves_material_unchanged(self): + p = MaterialDensity('Co', 8.9, 'Cobalt') + mw = p.molecular_weight.value + sld = p.sld.value + with self.assertRaises(ValueError): + p.chemical_structure = '###' + assert p.chemical_structure == 'Co' + assert_almost_equal(p.molecular_weight.value, mw) + assert_almost_equal(p.sld.value, sld) + + def test_sld_coupled_default_true(self): + p = MaterialDensity() + assert p.sld_coupled is True + assert p.sld.independent is False + assert p.isld.independent is False + + def test_decouple_keeps_values_and_detaches_density(self): + p = MaterialDensity(chemical_structure='Si', density=2.33) + coupled_sld = p.sld.value + p.sld_coupled = False + assert p.sld_coupled is False + assert p.sld.independent is True + assert p.isld.independent is True + assert_almost_equal(p.sld.value, coupled_sld) + # Density edits no longer propagate; sld is directly settable. + p.density.value = 9.99 + assert_almost_equal(p.sld.value, coupled_sld) + p.sld.value = 5.5 + assert p.sld.value == 5.5 + + def test_recouple_recomputes_and_discards_manual_sld(self): + p = MaterialDensity(chemical_structure='Si', density=2.33) + coupled_sld = p.sld.value + p.sld_coupled = False + p.sld.value = 5.5 + p.sld_coupled = True + assert p.sld_coupled is True + assert_almost_equal(p.sld.value, coupled_sld) + # And propagation is restored. + p.density.value = 4.66 + assert_almost_equal(p.sld.value, 2 * coupled_sld) + + def test_sld_coupled_setter_is_idempotent(self): + p = MaterialDensity() + p.sld_coupled = True # no-op, must not raise or rewire + p.sld_coupled = False + p.sld_coupled = False # no-op on the decoupled side too + assert p.sld_coupled is False + + def test_decoupled_sld_can_be_freed_for_fitting(self): + p = MaterialDensity() + p.sld_coupled = False + p.sld.fixed = False + free = p.get_fit_parameters() + assert any(parameter is p.sld for parameter in free) + + def test_dict_round_trip_decoupled_preserves_manual_sld(self): + p = MaterialDensity(chemical_structure='Si', density=2.33) + p.sld_coupled = False + p.sld.value = 7.25 + p.isld.value = -0.5 + p_dict = p.as_dict() + assert p_dict['sld_coupled'] is False + assert p_dict['sld'] == 7.25 + assert p_dict['isld'] == -0.5 + global_object.map._clear() + + q = MaterialDensity.from_dict(p_dict) + assert q.sld_coupled is False + assert_almost_equal(q.sld.value, 7.25) + assert_almost_equal(q.isld.value, -0.5) + # Still decoupled: density edits must not clobber the restored values. + q.density.value = 1.0 + assert_almost_equal(q.sld.value, 7.25) + + def test_dict_round_trip_coupled_carries_flag_but_no_sld(self): + p = MaterialDensity() + p_dict = p.as_dict() + assert p_dict['sld_coupled'] is True + assert 'sld' not in p_dict + assert 'isld' not in p_dict + + def test_from_dict_without_flag_restores_coupled(self): + """Project files predating the feature restore with current behavior.""" + p = MaterialDensity(chemical_structure='Si', density=2.33) + p_dict = {k: v for k, v in p.as_dict().items() if k != 'sld_coupled'} + global_object.map._clear() + + q = MaterialDensity.from_dict(p_dict) + assert q.sld_coupled is True + q.density.value = 4.66 + assert_almost_equal(q.sld.value, 2 * p.sld.value) + def test_density_mutation_propagates_after_round_trip(self): """Regression: after ``from_dict`` reattaches the saved ``_density`` Parameter, mutating it must propagate to ``sld`` / ``isld`` (which