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/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..3d7faccd --- /dev/null +++ b/docs/src/tutorials/simulation/bilayer.ipynb @@ -0,0 +1,696 @@ +{ + "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 solvent material — D₂O (heavy water) — which is used for all layers in the bilayer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76bd9056", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the solvent material\n", + "d2o = Material(sld=6.36, isld=0, name='D2O')" + ] + }, + { + "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", + "front_tail_layer = LayerAreaPerMolecule(\n", + " molecular_formula='C32D64',\n", + " thickness=16.0,\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", + " 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", + " 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", + ")\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: front_head, front_tail, back_tail, back_head\n", + "# All structural parameters are automatically constrained\n", + "print(f'Bilayer layers: {len(bilayer.layers)}')" + ] + }, + { + "cell_type": "markdown", + "id": "d9b60931", + "metadata": {}, + "source": [ + "### Verifying Constraints\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a607ca1d", + "metadata": {}, + "outputs": [], + "source": [ + "# 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} Ų')" + ] + }, + { + "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`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffb84f80", + "metadata": {}, + "outputs": [], + "source": [ + "# 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, + "id": "39802839", + "metadata": {}, + "outputs": [], + "source": [ + "# Conformal roughness: all layers share the same roughness value\n", + "# (controlled by the front head layer)\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 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 Substrate')\n", + "sio2_layer = Layer(material=sio2, thickness=15, roughness=3, name='SiO2')\n", + "\n", + "# D2O subphase (bulk water)\n", + "d2o_subphase = Layer(material=d2o, thickness=0, roughness=3, name='D2O Bulk')\n", + "\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(si_layer, name='Si Substrate'),\n", + " Multilayer(sio2_layer, name='SiO2'),\n", + " bilayer,\n", + " Multilayer(d2o_subphase, name='D2O Subphase'),\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": [ + "# 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\n", + "\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", + "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()\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()" + ] + }, + { + "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": [ + "# 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", + "# Get symmetric SLD profile\n", + "z_sym, sld_sym = model.interface().sld_profile(model.unique_name)\n", + "\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", + "# 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", + "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", + "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", + "\n", + "# Create H2O material for second contrast\n", + "h2o = Material(sld=-0.56, isld=0, name='H2O')\n", + "\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", + " 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 (same deuterated lipid, different solvent)\n", + "tail_layer_h2o = LayerAreaPerMolecule(\n", + " molecular_formula='C32D64',\n", + " thickness=16.0,\n", + " solvent=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.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('Reflectivity: Multiple Contrast Comparison')\n", + "plt.legend()\n", + "plt.grid(True, alpha=0.3)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "624b1d60", + "metadata": {}, + "source": [ + "## Summary\n", + "\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. **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 (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." + ] + } + ], + "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.12" + } + }, + "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 diff --git a/pyproject.toml b/pyproject.toml index 630422ca..0dde562c 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", @@ -35,7 +36,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", ] @@ -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/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..21c428f6 --- /dev/null +++ b/src/easyreflectometry/sample/assemblies/bilayer.py @@ -0,0 +1,491 @@ +from __future__ import annotations + +from typing import Any + +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 + +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. + + The bilayer structure is: Front Head - Front Tail - Back Tail - Back Head + + This assembly comes pre-populated with physically meaningful constraints: + - 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`_ + + .. _`bilayer documentation`: ../sample/assemblies_library.html#bilayer + """ + + def __init__( + self, + front_head_layer: LayerAreaPerMolecule | None = None, + front_tail_layer: LayerAreaPerMolecule | None = None, + back_head_layer: LayerAreaPerMolecule | None = None, + name: str = 'EasyBilayer', + unique_name: str | None = None, + constrain_heads: bool = True, + conformal_roughness: bool = True, + interface: Any = None, + ): + """Constructor. + + :param front_head_layer: Layer representing the front head part of the bilayer. + :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 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 layers if not provided + if front_head_layer is None: + front_head_layer = self._create_default_head_layer( + unique_name=unique_name, + name_suffix='Front', + interface=interface, + ) + + if front_tail_layer is None: + front_tail_layer = self._create_default_tail_layer( + unique_name=unique_name, + interface=interface, + ) + + # 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, + ) + + if back_head_layer is None: + back_head_layer = self._create_default_head_layer( + unique_name=unique_name, + name_suffix='Back', + interface=interface, + ) + + # Create layer collection: front_head, front_tail, back_tail, back_head + bilayer_layers = LayerCollection( + front_head_layer, + front_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 + + @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. + + 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( + dependency_expression='a', + dependency_map={'a': front_tail.thickness}, + ) + + # Constrain area per molecule + back_tail.area_per_molecule_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_tail.area_per_molecule_parameter}, + ) + + # Constrain solvent fraction + back_tail.solvent_fraction_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': front_tail.solvent_fraction_parameter}, + ) + + self._tail_constraints_setup = True + + @property + def front_head_layer(self) -> 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) -> LayerAreaPerMolecule: + """Get the front tail layer of the bilayer.""" + return self.layers[1] + + @property + def back_tail_layer(self) -> LayerAreaPerMolecule: + """Get the back tail layer of the bilayer.""" + return self.layers[2] + + @property + def back_head_layer(self) -> 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_parameter.make_dependent_on( + dependency_expression='a', + 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_parameter.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. + + 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. + :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_parameter.make_dependent_on( + dependency_expression='a', + 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_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.back_head_layer.area_per_molecule_parameter}, + ) + + if tail_area_per_molecule: + 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_parameter}, + ) + + if front_head_fraction: + self.front_head_layer.solvent_fraction_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_head_layer.solvent_fraction_parameter}, + ) + + if back_head_fraction: + self.back_head_layer.solvent_fraction_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.back_head_layer.solvent_fraction_parameter}, + ) + + if tail_fraction: + self.front_tail_layer.solvent_fraction_parameter.make_dependent_on( + dependency_expression='a', + dependency_map={'a': another_contrast.front_tail_layer.solvent_fraction_parameter}, + ) + + @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: list[str] | None = 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['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 + 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..b96e6a78 --- /dev/null +++ b/tests/sample/assemblies/test_bilayer.py @@ -0,0 +1,420 @@ +""" +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 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() + 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_matched_water = Material(sld=0, isld=0, name='Air Matched Water') + + 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_matched_water, + 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, + front_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 = 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 = 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_matched_water = Material(sld=0, isld=0, name='Air Matched Water') + + 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_matched_water, + 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, + front_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) 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