Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/python-ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ jobs:
name: unit-tests-job
flags: unittests
files: ./coverage-unit.xml
fail_ci_if_error: true
fail_ci_if_error: false
verbose: true
token: ${{ secrets.CODECOV_TOKEN }}
slug: EasyScience/EasyReflectometryLib
Expand Down
9 changes: 9 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,12 @@
# Version 1.6.0 (1 May 2026)

Add Mighell-based handling of non-positive-variance points in fitting (issue #256).
Non-positive-variance data points are no longer forcibly discarded; instead, a
hybrid objective applies a Mighell substitution for non-positive-variance points
while using standard weighted least squares for the rest. The previous masking
behavior is available via `objective='legacy_mask'`. New `objective` parameter on
`MultiFitter`, `fit()`, and `fit_single_data_set_1d()`.

# Version 1.3.3 (17 June 2025)

Added Chi^2 and fit status to fitting results.
Expand Down
9 changes: 9 additions & 0 deletions docs/src/api/api.rst
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,15 @@ Project provides a higher-level interface for managing models, experiments, and

project

Fitting
=======
Fitting helpers and objective functions.

.. toctree::
:maxdepth: 1

fitting

Assemblies
==========
Assemblies are collections of layers that are used to represent a specific physical setup.
Expand Down
112 changes: 112 additions & 0 deletions docs/src/api/fitting.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
Fitting
=======

.. currentmodule:: easyreflectometry.fitting

Objective functions and non-positive variance handling
-----------------------------------------------------

Check warning on line 7 in docs/src/api/fitting.rst

View workflow job for this annotation

GitHub Actions / build_documentation

Title underline too short.

:class:`MultiFitter` supports several objective modes for handling reflectometry
data during fitting, especially when measured variances are non-positive.

The default objective is ``hybrid``. This uses ordinary weighted least squares
for points with positive variance and applies a Mighell-style substitution only
to points whose variance is non-positive. The older ``legacy_mask`` mode
drops non-positive-variance points before fitting. The ``mighell`` mode applies the
Mighell transform to every point.

Mighell objective
~~~~~~~~~~~~~~~~~

The full ``mighell`` objective follows the algebraic form of the
``chi^2_gamma`` statistic described by Mighell for Poisson-distributed count
data:

.. math::

\chi^2_\gamma =
\sum_i \frac{[n_i + \min(n_i, 1) - m_i]^2}{n_i + 1}

where ``n_i`` are observed counts and ``m_i`` are model values.

In EasyReflectometry this is implemented as a weighted least-squares problem.
For each observed value ``y_i`` the fitted target is shifted to

.. math::

y_{\mathrm{eff},i} = y_i + \min(y_i, 1)

and the effective uncertainty is

.. math::

\sigma_i = \sqrt{y_i + 1}

so the minimized objective is

.. math::

\sum_i \left(\frac{y_{\mathrm{eff},i} - f_i}{\sigma_i}\right)^2 =
\sum_i \frac{[y_i + \min(y_i, 1) - f_i]^2}{y_i + 1}

This is the same algebraic form as Mighell's statistic, with the model value
``f_i`` replacing ``m_i``.

Scope and interpretation
~~~~~~~~~~~~~~~~~~~~~~~~

Mighell's statistic was derived for Poisson-distributed count data. In
reflectometry workflows, the fitted values are usually normalized
reflectivities or intensities rather than raw counts. They may already have
been processed, scaled, background-corrected, or otherwise transformed before
they reach the fitter.

This distinction matters when interpreting the result. The full ``mighell``
objective is not only a reweighting of residuals; it also changes the fitted
target from ``y`` to ``y + min(y, 1)``. For values between zero and one, this
can substantially increase the target value. A fit can therefore have a good
Mighell objective value while looking poorer against the originally plotted
reflectivity curve, or while having a worse classical chi-square.

For reflectometry data, ``hybrid`` is generally the recommended compromise:
it preserves ordinary weighted least-squares behavior where positive variances are
available, while still allowing non-positive-variance points to contribute through the
Mighell-style substitution.

Objective modes
~~~~~~~~~~~~~~~

``hybrid``
Default. Use standard weighted least squares for points with positive
variance and apply the Mighell substitution only where variance is
non-positive.

``mighell``
Apply the Mighell transform to all points. The reported objective chi-square
is evaluated in transformed objective space and should not be interpreted as
a classical chi-square against the original reflectivity values.

``legacy_mask``
Remove non-positive-variance points before fitting and use standard weighted least
squares for the remaining points.

``auto``
Alias for ``hybrid``.

Fit metrics
~~~~~~~~~~~

The fitter exposes both objective-space and classical fit metrics after fitting.
``objective_chi2`` and ``objective_reduced_chi`` describe the minimized
objective, which may include transformed targets under ``hybrid`` or
``mighell``. ``classical_chi2`` and ``classical_reduced_chi`` are computed
against the original observed reflectivity values using only points with
positive variance.

API reference
-------------

.. automodule:: easyreflectometry.fitting
:members:
:undoc-members:
:show-inheritance:
117 changes: 111 additions & 6 deletions docs/src/tutorials/fitting/simple_fitting.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,9 @@
"source": [
"%matplotlib inline\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import pooch\n",
"import refl1d\n",
"import refnx\n",
"\n",
"import easyreflectometry\n",
Expand Down Expand Up @@ -63,7 +65,8 @@
"outputs": [],
"source": [
"print(f'easyreflectometry: {easyreflectometry.__version__}')\n",
"print(f'refnx: {refnx.__version__}')"
"print(f'refnx: {refnx.__version__}')\n",
"print(f'refl1d: {refl1d.__version__}')"
]
},
{
Expand Down Expand Up @@ -395,7 +398,7 @@
"## Choosing our calculation engine\n",
"\n",
"The `easyreflectometry` package enables the calculation of the reflectometry profile using either [*refnx*](https://refnx.readthedocs.io/) or [*Refl1D*](https://refl1d.readthedocs.io/en/latest/).\n",
"For this tutorial, we will stick to the current default, which is *refnx*. \n",
"We will first run the fit with the current default, *refnx*, and then rebuild the same starting model with *Refl1D* to compare the results.\n",
"The calculator must be created and associated with the model that we are to fit. "
]
},
Expand Down Expand Up @@ -464,7 +467,8 @@
"metadata": {},
"outputs": [],
"source": [
"analysed = fitter.fit(data)"
"analysed_refnx = fitter.fit(data)\n",
"analysed = analysed_refnx"
]
},
{
Expand Down Expand Up @@ -513,18 +517,119 @@
"model"
]
},
{
"cell_type": "markdown",
"id": "2b2fe558",
"metadata": {},
"source": [
"## Repeating the optimisation with Refl1D\n",
"\n",
"To compare backends fairly, we rebuild the same initial model and then switch the calculator to `refl1d`.\n",
"This ensures that the second optimisation starts from the same parameter guesses rather than from the already-optimised `refnx` result."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "341abc6f",
"metadata": {},
"outputs": [],
"source": [
"print(f'refl1d: {refl1d.__version__}')\n",
"\n",
"si_refl1d = Material(sld=2.07, isld=0, name='Si')\n",
"sio2_refl1d = Material(sld=3.47, isld=0, name='SiO2')\n",
"film_refl1d = Material(sld=2.0, isld=0, name='Film')\n",
"d2o_refl1d = Material(sld=6.36, isld=0, name='D2O')\n",
"\n",
"si_layer_refl1d = Layer(material=si_refl1d, thickness=0, roughness=0, name='Si layer')\n",
"sio2_layer_refl1d = Layer(material=sio2_refl1d, thickness=30, roughness=3, name='SiO2 layer')\n",
"film_layer_refl1d = Layer(material=film_refl1d, thickness=250, roughness=3, name='Film Layer')\n",
"subphase_refl1d = Layer(material=d2o_refl1d, thickness=0, roughness=3, name='D2O Subphase')\n",
"\n",
"superphase_refl1d = Multilayer([si_layer_refl1d, sio2_layer_refl1d], name='Si/SiO2 Superphase')\n",
"sample_refl1d = Sample(superphase_refl1d, Multilayer(film_layer_refl1d), Multilayer(subphase_refl1d), name='Film Structure')\n",
"\n",
"resolution_function_refl1d = PercentageFwhm(0.02)\n",
"model_refl1d = Model(\n",
" sample=sample_refl1d,\n",
" scale=1,\n",
" background=1e-6,\n",
" resolution_function=resolution_function_refl1d,\n",
" name='Film Model (Refl1D)'\n",
")\n",
"\n",
"sio2_layer_refl1d.thickness.bounds = (15, 50)\n",
"film_layer_refl1d.thickness.bounds = (200, 300)\n",
"sio2_layer_refl1d.roughness.bounds = (1, 15)\n",
"film_layer_refl1d.roughness.bounds = (1, 15)\n",
"subphase_refl1d.roughness.bounds = (1, 15)\n",
"film_layer_refl1d.material.sld.bounds = (0.1, 3)\n",
"model_refl1d.background.bounds = (1e-8, 1e-5)\n",
"model_refl1d.scale.bounds = (0.5, 1.5)\n",
"\n",
"interface_refl1d = CalculatorFactory()\n",
"interface_refl1d.switch('refl1d')\n",
"model_refl1d.interface = interface_refl1d\n",
"\n",
"print(interface_refl1d.current_interface.name)"
]
},
{
"cell_type": "markdown",
"id": "2f764d79",
"metadata": {},
"source": [
"We can now fit the same dataset with `refl1d` and compare the fitted curve and reduced chi-squared with the earlier `refnx` result."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70229cf9",
"metadata": {},
"outputs": [],
"source": [
"fitter_refl1d = MultiFitter(model_refl1d)\n",
"analysed_refl1d = fitter_refl1d.fit(data)\n",
"\n",
"print(f'refnx reduced chi: {analysed_refnx[\"reduced_chi\"]:.4f}')\n",
"print(f'refl1d reduced chi: {analysed_refl1d[\"reduced_chi\"]:.4f}')\n",
"\n",
"qz = data['coords']['Qz_0'].values\n",
"reflectivity = data['data']['R_0'].values\n",
"uncertainty = data['data']['R_0'].variances**0.5\n",
"\n",
"plt.figure(figsize=(8, 5))\n",
"plt.errorbar(qz, reflectivity, yerr=uncertainty, fmt='o', markersize=3, color='black', alpha=0.6, label='Data')\n",
"plt.plot(qz, analysed_refnx['R_0_model'].values, linewidth=2.5, label='refnx fit')\n",
"plt.plot(qz, analysed_refl1d['R_0_model'].values, linewidth=2.5, label='refl1d fit')\n",
"plt.yscale('log')\n",
"plt.xlabel(r'$Q_z$ (1/Å)')\n",
"plt.ylabel('Reflectivity')\n",
"plt.grid(True, which='both', alpha=0.25)\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"print('refnx model')\n",
"print(model)\n",
"print()\n",
"print('refl1d model')\n",
"print(model_refl1d)"
]
},
{
"cell_type": "markdown",
"id": "df6db8c3-6515-478b-bd2e-aac967579231",
"metadata": {},
"source": [
"We note here that the results obtained are very similar to those from the [*refnx* tutorial](https://refnx.readthedocs.io/en/latest/getting_started.html#Fitting-a-neutron-reflectometry-dataset), which is hardly surprising given that we have used the *refnx* engine in this example."
"We note here that the results obtained with [*refnx*](https://refnx.readthedocs.io/en/latest/) and [*Refl1D*](https://refl1d.readthedocs.io/en/latest/) are very similar for this simple slab model, with only small backend-dependent differences in the fitted curve and optimised parameters."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "easyref",
"display_name": ".venv (3.11.9)",
"language": "python",
"name": "python3"
},
Expand All @@ -538,7 +643,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.9"
"version": "3.11.9"
}
},
"nbformat": 4,
Expand Down
Loading
Loading