From 3b825d7bae5b80315905f26e261e175a6d2bede5 Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:50:45 +0800 Subject: [PATCH 01/75] chore: establish package structure and data policy --- .gitignore | 15 ++++++++++++++- README.md | 32 ++++++++++++++++++++++++++++++++ docs/data_policy.md | 11 +++++++++++ pyproject.toml | 28 ++++++++++++++++++++++++++-- src/xps_fitting/__init__.py | 3 +++ tests/test_package.py | 4 ++++ 6 files changed, 90 insertions(+), 3 deletions(-) create mode 100644 README.md create mode 100644 docs/data_policy.md create mode 100644 src/xps_fitting/__init__.py create mode 100644 tests/test_package.py diff --git a/.gitignore b/.gitignore index 35c3b19..acbfa47 100644 --- a/.gitignore +++ b/.gitignore @@ -16,4 +16,17 @@ __pycache__/ .ipynb_checkpoints/ # macOS -.DS_Store \ No newline at end of file +.DS_Store + +# Experimental and generated data (small curated fixtures remain allowed) +*.vgd +*.VGD +raw_data/ +processed_data/ +processed_xps/ +outputs/ +figures/ +reports/ +build/ +dist/ +*.egg-info/ diff --git a/README.md b/README.md new file mode 100644 index 0000000..0d4e974 --- /dev/null +++ b/README.md @@ -0,0 +1,32 @@ +# XPS fitting workbench + +This project is a reproducible, chemistry-constrained X-ray photoelectron +spectroscopy fitting engine. It optimises user-proposed chemical components and +supports comparison of alternative hypotheses; it does not discover chemical +assignments or endorse a model merely because it has more peaks. + +> Peak fitting requires chemical judgement. Fit statistics alone cannot establish +> a chemically correct assignment. + +## Status + +Phase 1 is under active development. It establishes numerical models, diagnostics, +exports, and a stable result interface. Colourful publication figures are deferred +to Phase 2; Phase 1 plotting is diagnostic only. + +## Installation + +```bash +python -m venv .venv +source .venv/bin/activate +python -m pip install -e '.[dev]' +pytest +``` + +The project uses a `src/xps_fitting/` package, `tests/`, `examples/`, `configs/`, +and `docs/`. Raw experimental data, generated outputs, reports, and figures are +ignored. Deliberately curated small text fixtures may be committed under +`tests/data/` or `examples/data/`; see [the data policy](docs/data_policy.md). + +The legacy `xps_vgd_utils.py` workflow remains in place while package I/O support +is developed. diff --git a/docs/data_policy.md b/docs/data_policy.md new file mode 100644 index 0000000..b597e9f --- /dev/null +++ b/docs/data_policy.md @@ -0,0 +1,11 @@ +# Data policy + +Only small, anonymised, scientifically useful text fixtures belong in `tests/data/` +or `examples/data/`. A reduced experimental fixture must remove private paths and +sample identifiers and document the source range, columns, downsampling, and any +normalisation. Raw VGD files and bulk `raw_data/` or `processed_data/` directories +are excluded because they may contain private metadata and are unsuitable for Git. + +The raw files already tracked in this repository predate this policy. They are left +untouched, are not copied into tests, and should be reviewed by the data owner. +Deterministic synthetic spectra are the default for tests and examples. diff --git a/pyproject.toml b/pyproject.toml index 2224004..d902a93 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,29 @@ [project] name = "xps-fitting-workbench" version = "0.1.0" -requires-python = ">=3.14" -dependencies = [] +description = "Reproducible, chemistry-constrained XPS peak fitting" +readme = "README.md" +requires-python = ">=3.10" +dependencies = [ + "lmfit>=1.2", + "matplotlib>=3.7", + "numpy>=1.24", + "openpyxl>=3.0", + "pandas>=2.0", + "scipy>=1.10", +] + +[project.optional-dependencies] +dev = ["pytest>=8"] + +[build-system] +requires = ["setuptools>=68"] +build-backend = "setuptools.build_meta" + +[tool.setuptools.packages.find] +where = ["src"] + +[tool.pytest.ini_options] +addopts = "-ra" +testpaths = ["tests"] +pythonpath = ["src", "."] diff --git a/src/xps_fitting/__init__.py b/src/xps_fitting/__init__.py new file mode 100644 index 0000000..42b3897 --- /dev/null +++ b/src/xps_fitting/__init__.py @@ -0,0 +1,3 @@ +"""Chemistry-constrained XPS fitting tools.""" + +__version__ = "0.1.0" diff --git a/tests/test_package.py b/tests/test_package.py new file mode 100644 index 0000000..394aa06 --- /dev/null +++ b/tests/test_package.py @@ -0,0 +1,4 @@ +def test_package_imports() -> None: + import xps_fitting + + assert xps_fitting.__version__ == "0.1.0" From 35bdfd4954da29ef6f81387d5cc341369ee68d12 Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:51:47 +0800 Subject: [PATCH 02/75] feat: add validated spectrum input and example data workflow --- README.md | 16 +++++++++++- examples/load_spectrum.py | 9 +++++++ src/xps_fitting/__init__.py | 4 +++ src/xps_fitting/io.py | 45 +++++++++++++++++++++++++++++++++ src/xps_fitting/io_vgd.py | 25 +++++++++++++++++++ src/xps_fitting/spectrum.py | 50 +++++++++++++++++++++++++++++++++++++ tests/test_io.py | 22 ++++++++++++++++ 7 files changed, 170 insertions(+), 1 deletion(-) create mode 100644 examples/load_spectrum.py create mode 100644 src/xps_fitting/io.py create mode 100644 src/xps_fitting/io_vgd.py create mode 100644 src/xps_fitting/spectrum.py create mode 100644 tests/test_io.py diff --git a/README.md b/README.md index 0d4e974..b5b5348 100644 --- a/README.md +++ b/README.md @@ -29,4 +29,18 @@ ignored. Deliberately curated small text fixtures may be committed under `tests/data/` or `examples/data/`; see [the data policy](docs/data_policy.md). The legacy `xps_vgd_utils.py` workflow remains in place while package I/O support -is developed. +is developed. CSV and XLSX tables, pandas DataFrames, and VGD through the optional +legacy adapter are supported. All input is cleaned, duplicate energies are averaged, +and fitting arrays are ascending regardless of acquisition order; the original order +is recorded in metadata. + +```python +from xps_fitting.io import read_csv +spectrum = read_csv("spectrum.csv", region="C 1s", sample_name="sample") +``` + +Existing notebooks may continue importing `xps_vgd_utils`. New code can use +`xps_fitting.io_vgd.read_vgd`; it requires the separately installed `vgd-reader` +package. Examples and tests use deterministic synthetic data rather than copying +the tracked experimental files. Run `python examples/load_spectrum.py` after an +editable install. diff --git a/examples/load_spectrum.py b/examples/load_spectrum.py new file mode 100644 index 0000000..8450067 --- /dev/null +++ b/examples/load_spectrum.py @@ -0,0 +1,9 @@ +"""Load a deliberately messy synthetic table through the public I/O API.""" + +import pandas as pd + +from xps_fitting.io import spectrum_from_dataframe + +frame = pd.DataFrame({"energy": [286.0, 285.0, 285.0, "bad", 284.0], "counts": [1, 3, 5, 9, 2]}) +spectrum = spectrum_from_dataframe(frame, "energy", "counts", region="C 1s", sample_name="synthetic") +print(f"Loaded {spectrum.binding_energy.size} ordered points: {spectrum.binding_energy}") diff --git a/src/xps_fitting/__init__.py b/src/xps_fitting/__init__.py index 42b3897..e68bd0e 100644 --- a/src/xps_fitting/__init__.py +++ b/src/xps_fitting/__init__.py @@ -1,3 +1,7 @@ """Chemistry-constrained XPS fitting tools.""" +from .spectrum import Spectrum + __version__ = "0.1.0" + +__all__ = ["Spectrum", "__version__"] diff --git a/src/xps_fitting/io.py b/src/xps_fitting/io.py new file mode 100644 index 0000000..22218e2 --- /dev/null +++ b/src/xps_fitting/io.py @@ -0,0 +1,45 @@ +"""Tabular spectrum input.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd + +from .spectrum import Spectrum + + +def spectrum_from_dataframe( + frame: pd.DataFrame, + energy_column: str = "binding_energy_eV", + intensity_column: str = "intensity", + *, region: str = "", sample_name: str = "", source_file: str | None = None, + metadata: dict[str, Any] | None = None, +) -> Spectrum: + """Clean a table, average duplicate energies, and return ascending arrays.""" + if energy_column not in frame or intensity_column not in frame: + raise ValueError(f"required columns: {energy_column!r}, {intensity_column!r}") + original_order = "ascending" if frame[energy_column].is_monotonic_increasing else ( + "descending" if frame[energy_column].is_monotonic_decreasing else "unordered" + ) + table = pd.DataFrame({ + "energy": pd.to_numeric(frame[energy_column], errors="coerce"), + "intensity": pd.to_numeric(frame[intensity_column], errors="coerce"), + }).dropna() + table = table[np.isfinite(table["energy"]) & np.isfinite(table["intensity"])] + table = table.groupby("energy", as_index=False, sort=True)["intensity"].mean() + info = dict(metadata or {}) + info.update({"original_order": original_order, "input_rows": len(frame), "clean_rows": len(table)}) + return Spectrum(table.energy.to_numpy(), table.intensity.to_numpy(), region, sample_name, source_file, info) + + +def read_csv(path: str | Path, **kwargs: Any) -> Spectrum: + path = Path(path) + return spectrum_from_dataframe(pd.read_csv(path), source_file=str(path), **kwargs) + + +def read_xlsx(path: str | Path, *, sheet_name: str | int = 0, **kwargs: Any) -> Spectrum: + path = Path(path) + return spectrum_from_dataframe(pd.read_excel(path, sheet_name=sheet_name), source_file=str(path), **kwargs) diff --git a/src/xps_fitting/io_vgd.py b/src/xps_fitting/io_vgd.py new file mode 100644 index 0000000..db35275 --- /dev/null +++ b/src/xps_fitting/io_vgd.py @@ -0,0 +1,25 @@ +"""Optional adapter for the repository's existing VGD reader workflow.""" + +from __future__ import annotations + +from pathlib import Path + +from .io import spectrum_from_dataframe +from .spectrum import Spectrum + + +def read_vgd(path: str | Path, *, spectrum_index: int = 0, intensity_column: str = "corrected_intensity") -> Spectrum: + try: + from xps_vgd_utils import read_vgd_to_dataframe + except ImportError as exc: + raise ImportError("VGD input requires the optional vgd-reader package and xps_vgd_utils.py") from exc + frame = read_vgd_to_dataframe(path) + selected = frame[frame["spectrum_index"] == spectrum_index] + if selected.empty: + raise IndexError(f"spectrum_index {spectrum_index} is unavailable") + row = selected.iloc[0] + return spectrum_from_dataframe( + selected, "binding_energy_eV", intensity_column, + region=str(row.get("core_level", "")), sample_name=str(row.get("sample_id", "")), source_file=str(path), + metadata={"vgd_spectrum_index": spectrum_index}, + ) diff --git a/src/xps_fitting/spectrum.py b/src/xps_fitting/spectrum.py new file mode 100644 index 0000000..e06e36e --- /dev/null +++ b/src/xps_fitting/spectrum.py @@ -0,0 +1,50 @@ +"""Validated, numerically safe spectrum representation.""" + +from __future__ import annotations + +from dataclasses import dataclass, field, replace +from pathlib import Path +from typing import Any + +import numpy as np + + +@dataclass(frozen=True) +class Spectrum: + binding_energy: np.ndarray + intensity: np.ndarray + region: str = "" + sample_name: str = "" + source_file: str | None = None + metadata: dict[str, Any] = field(default_factory=dict) + normalised_intensity: np.ndarray | None = None + acquisition_metadata: dict[str, Any] = field(default_factory=dict) + + def __post_init__(self) -> None: + energy = np.asarray(self.binding_energy, dtype=float) + intensity = np.asarray(self.intensity, dtype=float) + if energy.ndim != 1 or intensity.ndim != 1 or energy.size != intensity.size: + raise ValueError("binding_energy and intensity must be equal-length 1D arrays") + if energy.size < 2 or not np.all(np.isfinite(energy)) or not np.all(np.isfinite(intensity)): + raise ValueError("a spectrum needs at least two finite points") + if np.any(np.diff(energy) <= 0): + raise ValueError("binding_energy must be strictly ascending for fitting") + object.__setattr__(self, "binding_energy", energy) + object.__setattr__(self, "intensity", intensity) + if self.normalised_intensity is not None: + norm = np.asarray(self.normalised_intensity, dtype=float) + if norm.shape != energy.shape: + raise ValueError("normalised_intensity must match binding_energy") + object.__setattr__(self, "normalised_intensity", norm) + + def crop(self, minimum: float, maximum: float) -> "Spectrum": + lo, hi = sorted((minimum, maximum)) + mask = (self.binding_energy >= lo) & (self.binding_energy <= hi) + if mask.sum() < 2: + raise ValueError("crop interval contains fewer than two points") + normalised = None if self.normalised_intensity is None else self.normalised_intensity[mask] + return replace(self, binding_energy=self.binding_energy[mask], intensity=self.intensity[mask], normalised_intensity=normalised) + + @property + def source_path(self) -> Path | None: + return None if self.source_file is None else Path(self.source_file) diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..e83d1bb --- /dev/null +++ b/tests/test_io.py @@ -0,0 +1,22 @@ +import numpy as np +import pandas as pd + +from xps_fitting.io import read_csv, read_xlsx, spectrum_from_dataframe + + +def test_dataframe_cleaning_order_duplicates_and_crop() -> None: + frame = pd.DataFrame({"E": [3, 2, 2, "bad", 1, np.nan], "I": [4, 2, 4, 9, 1, 5]}) + spectrum = spectrum_from_dataframe(frame, "E", "I") + np.testing.assert_allclose(spectrum.binding_energy, [1, 2, 3]) + np.testing.assert_allclose(spectrum.intensity, [1, 3, 4]) + assert spectrum.metadata["original_order"] == "unordered" + np.testing.assert_allclose(spectrum.crop(1, 2).binding_energy, [1, 2]) + + +def test_csv_and_xlsx(tmp_path) -> None: + frame = pd.DataFrame({"binding_energy_eV": [2, 1], "intensity": [20, 10]}) + csv_path, xlsx_path = tmp_path / "s.csv", tmp_path / "s.xlsx" + frame.to_csv(csv_path, index=False) + frame.to_excel(xlsx_path, index=False) + assert read_csv(csv_path).metadata["original_order"] == "descending" + np.testing.assert_allclose(read_xlsx(xlsx_path).intensity, [10, 20]) From dae026297097e7639264a98c081dea07adc00c87 Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:52:57 +0800 Subject: [PATCH 03/75] feat: implement XPS backgrounds line shapes and constraints --- README.md | 18 ++++++++++++ src/xps_fitting/backgrounds.py | 30 +++++++++++++++++++ src/xps_fitting/configuration.py | 50 ++++++++++++++++++++++++++++++++ src/xps_fitting/constraints.py | 30 +++++++++++++++++++ src/xps_fitting/lineshapes.py | 31 ++++++++++++++++++++ tests/test_backgrounds.py | 13 +++++++++ tests/test_constraints.py | 14 +++++++++ tests/test_lineshapes.py | 17 +++++++++++ 8 files changed, 203 insertions(+) create mode 100644 src/xps_fitting/backgrounds.py create mode 100644 src/xps_fitting/configuration.py create mode 100644 src/xps_fitting/constraints.py create mode 100644 src/xps_fitting/lineshapes.py create mode 100644 tests/test_backgrounds.py create mode 100644 tests/test_constraints.py create mode 100644 tests/test_lineshapes.py diff --git a/README.md b/README.md index b5b5348..b931ea1 100644 --- a/README.md +++ b/README.md @@ -44,3 +44,21 @@ Existing notebooks may continue importing `xps_vgd_utils`. New code can use package. Examples and tests use deterministic synthetic data rather than copying the tracked experimental files. Run `python examples/load_spectrum.py` after an editable install. + +## Physical model + +Implemented backgrounds are linear and iterative Shirley. Tougaard is deferred +until a responsible convention and validation dataset are selected. Gaussian, +Lorentzian, pseudo-Voigt, and true Voigt functions are area-normalised; widths are +FWHM (the true Voigt accepts separate Gaussian and Lorentzian FWHM values). + +Peak configuration supports non-negative/bounded areas, fixed or bounded values, +shared width and fraction groups, fixed centre offsets, fixed area ratios, and +spin-orbit doublets. For example, `cl2p_doublet("Cl", 200.0, 1000.0)` creates a +configurable 1.6 eV-separated pair with 2:1 area, shared FWHM, and shared mixing. +Bounds and assignments remain the analyst's responsibility. + +KherveFitting source was not found locally. Its named GL/LA functions can use +product or asymmetric conventions, whereas lmfit's `PseudoVoigtModel` is a linear +mixture of area-normalised Gaussian and Lorentzian profiles with a shared FWHM. +No numerical equivalence is claimed and no Kherve-compatible peak shape is exposed. diff --git a/src/xps_fitting/backgrounds.py b/src/xps_fitting/backgrounds.py new file mode 100644 index 0000000..41e8ac3 --- /dev/null +++ b/src/xps_fitting/backgrounds.py @@ -0,0 +1,30 @@ +"""XPS background estimators.""" + +from __future__ import annotations + +import numpy as np + + +def linear(x: np.ndarray, y_start: float, y_end: float) -> np.ndarray: + x = np.asarray(x, dtype=float) + if x[-1] == x[0]: + raise ValueError("energy range must be nonzero") + return y_start + (y_end - y_start) * (x - x[0]) / (x[-1] - x[0]) + + +def shirley(intensity: np.ndarray, *, max_iter: int = 200, tolerance: float = 1e-7) -> np.ndarray: + """Iterative discrete Shirley background for data in either acquisition order.""" + y = np.asarray(intensity, dtype=float) + if y.ndim != 1 or y.size < 2 or not np.all(np.isfinite(y)): + raise ValueError("intensity must be a finite 1D array") + background = np.linspace(y[0], y[-1], y.size) + for _ in range(max_iter): + positive = np.maximum(y - background, 0) + cumulative = np.cumsum(positive[::-1])[::-1] + if cumulative[0] <= np.finfo(float).eps: + return np.linspace(y[0], y[-1], y.size) + updated = y[-1] + (y[0] - y[-1]) * cumulative / cumulative[0] + if np.max(np.abs(updated - background)) <= tolerance * max(1.0, np.ptp(y)): + return updated + background = updated + return background diff --git a/src/xps_fitting/configuration.py b/src/xps_fitting/configuration.py new file mode 100644 index 0000000..c746fc1 --- /dev/null +++ b/src/xps_fitting/configuration.py @@ -0,0 +1,50 @@ +"""Validated fitting configuration dataclasses.""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass, field +from typing import Any + + +@dataclass +class PeakConfig: + label: str + centre: float + centre_bounds: tuple[float, float] + area: float + area_bounds: tuple[float, float] = (0.0, float("inf")) + fwhm: float = 1.4 + fwhm_bounds: tuple[float, float] = (0.2, 3.0) + line_shape: str = "pseudo_voigt" + fraction: float = 0.5 + fraction_bounds: tuple[float, float] = (0.0, 1.0) + fixed: tuple[str, ...] = () + width_group: str | None = None + fraction_group: str | None = None + centre_offset_from: tuple[str, float] | None = None + area_ratio_to: tuple[str, float] | None = None + + def __post_init__(self) -> None: + for name, value, bounds in (("centre", self.centre, self.centre_bounds), ("area", self.area, self.area_bounds), ("fwhm", self.fwhm, self.fwhm_bounds), ("fraction", self.fraction, self.fraction_bounds)): + if bounds[0] > value or value > bounds[1]: + raise ValueError(f"{self.label} {name} is outside its bounds") + if self.area < 0 or self.fwhm <= 0: + raise ValueError("area must be non-negative and FWHM positive") + + +@dataclass +class FitConfig: + name: str + region: str + peaks: list[PeakConfig] + background: str = "linear" + robust_loss: str = "linear" + multistart: int = 1 + random_seed: int = 0 + release_fraction: bool = False + max_background_iterations: int = 3 + width_penalty: float = 0.0 + metadata: dict[str, Any] = field(default_factory=dict) + + def to_dict(self) -> dict[str, Any]: + return asdict(self) diff --git a/src/xps_fitting/constraints.py b/src/xps_fitting/constraints.py new file mode 100644 index 0000000..110f356 --- /dev/null +++ b/src/xps_fitting/constraints.py @@ -0,0 +1,30 @@ +"""Helpers for common chemistry/physics parameter relationships.""" + +from __future__ import annotations + +from .configuration import PeakConfig + + +def cl2p_doublet( + label: str, centre_32: float, area_32: float, *, separation: float = 1.6, + fwhm: float = 1.2, fraction: float = 0.5, +) -> list[PeakConfig]: + """Create a 2p3/2:2p1/2 doublet with shared shape/width and a 2:1 area ratio.""" + group = f"{label}_doublet" + main = PeakConfig(f"{label}_2p3/2", centre_32, (centre_32 - 0.5, centre_32 + 0.5), area_32, fwhm=fwhm, width_group=group, fraction=fraction, fraction_group=group) + partner = PeakConfig( + f"{label}_2p1/2", centre_32 + separation, (centre_32 + separation - 0.5, centre_32 + separation + 0.5), area_32 / 2, + fwhm=fwhm, width_group=group, fraction=fraction, fraction_group=group, + centre_offset_from=(main.label, separation), area_ratio_to=(main.label, 0.5), + ) + return [main, partner] + + +def validate_links(peaks: list[PeakConfig]) -> None: + labels = {peak.label for peak in peaks} + if len(labels) != len(peaks): + raise ValueError("peak labels must be unique") + for peak in peaks: + for link in (peak.centre_offset_from, peak.area_ratio_to): + if link and link[0] not in labels: + raise ValueError(f"{peak.label} references unknown peak {link[0]}") diff --git a/src/xps_fitting/lineshapes.py b/src/xps_fitting/lineshapes.py new file mode 100644 index 0000000..3887219 --- /dev/null +++ b/src/xps_fitting/lineshapes.py @@ -0,0 +1,31 @@ +"""Area-normalised peak line shapes using FWHM throughout.""" + +from __future__ import annotations + +import numpy as np +from scipy.special import voigt_profile + + +def gaussian(x: np.ndarray, area: float, centre: float, fwhm: float) -> np.ndarray: + sigma = fwhm / (2 * np.sqrt(2 * np.log(2))) + return area * np.exp(-0.5 * ((x - centre) / sigma) ** 2) / (sigma * np.sqrt(2 * np.pi)) + + +def lorentzian(x: np.ndarray, area: float, centre: float, fwhm: float) -> np.ndarray: + gamma = fwhm / 2 + return area * gamma / (np.pi * ((x - centre) ** 2 + gamma**2)) + + +def pseudo_voigt(x: np.ndarray, area: float, centre: float, fwhm: float, fraction: float = 0.5) -> np.ndarray: + if not 0 <= fraction <= 1: + raise ValueError("fraction must lie in [0, 1]") + return fraction * lorentzian(x, area, centre, fwhm) + (1 - fraction) * gaussian(x, area, centre, fwhm) + + +def voigt(x: np.ndarray, area: float, centre: float, fwhm_g: float, fwhm_l: float) -> np.ndarray: + sigma = fwhm_g / (2 * np.sqrt(2 * np.log(2))) + gamma = fwhm_l / 2 + return area * voigt_profile(x - centre, sigma, gamma) + + +LINE_SHAPES = {"gaussian": gaussian, "lorentzian": lorentzian, "pseudo_voigt": pseudo_voigt, "voigt": voigt} diff --git a/tests/test_backgrounds.py b/tests/test_backgrounds.py new file mode 100644 index 0000000..d761659 --- /dev/null +++ b/tests/test_backgrounds.py @@ -0,0 +1,13 @@ +import numpy as np + +from xps_fitting.backgrounds import linear, shirley + + +def test_linear_endpoints_and_shirley_stability() -> None: + x = np.linspace(0, 10, 101) + np.testing.assert_allclose(linear(x, 2, 5)[[0, -1]], [2, 5]) + y = 2 + np.exp(-((x - 5) / 1.2) ** 2) + first = shirley(y) + np.testing.assert_allclose(first, shirley(y)) + assert np.all(np.isfinite(first)) + assert first[0] == y[0] and first[-1] == y[-1] diff --git a/tests/test_constraints.py b/tests/test_constraints.py new file mode 100644 index 0000000..85918c6 --- /dev/null +++ b/tests/test_constraints.py @@ -0,0 +1,14 @@ +import pytest + +from xps_fitting.configuration import PeakConfig +from xps_fitting.constraints import cl2p_doublet, validate_links + + +def test_nonnegative_area_and_valid_links() -> None: + with pytest.raises(ValueError): + PeakConfig("bad", 1, (0, 2), -1) + peaks = cl2p_doublet("Cl", 200.0, 100.0) + validate_links(peaks) + assert peaks[1].centre == pytest.approx(peaks[0].centre + 1.6) + assert peaks[1].area == pytest.approx(peaks[0].area / 2) + assert peaks[0].width_group == peaks[1].width_group diff --git a/tests/test_lineshapes.py b/tests/test_lineshapes.py new file mode 100644 index 0000000..adab69c --- /dev/null +++ b/tests/test_lineshapes.py @@ -0,0 +1,17 @@ +import numpy as np + +from xps_fitting.lineshapes import gaussian, lorentzian, pseudo_voigt, voigt + + +def test_requested_area_is_recovered() -> None: + x = np.linspace(-100, 100, 200_001) + for curve in (gaussian(x, 7, 0, 2), lorentzian(x, 7, 0, 2), pseudo_voigt(x, 7, 0, 2, 0.3), voigt(x, 7, 0, 2, 1)): + assert np.trapz(curve, x) == pytest.approx(7, rel=0.01) + + +def test_gaussian_and_lorentzian_fwhm() -> None: + for function in (gaussian, lorentzian): + assert function(np.array([1.0]), 1, 0, 2)[0] == pytest.approx(function(np.array([0.0]), 1, 0, 2)[0] / 2) + + +import pytest From a203c3289f0265b9155c54455ccd58e20f635084 Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:54:57 +0800 Subject: [PATCH 04/75] feat: add constrained optimiser diagnostics and FitResult API --- README.md | 20 +++++ src/xps_fitting/__init__.py | 3 +- src/xps_fitting/diagnostics.py | 19 +++++ src/xps_fitting/models.py | 41 +++++++++++ src/xps_fitting/optimiser.py | 118 ++++++++++++++++++++++++++++++ src/xps_fitting/plotting.py | 26 +++++++ src/xps_fitting/result.py | 42 +++++++++++ tests/test_doublet_constraints.py | 20 +++++ tests/test_synthetic_fit.py | 24 ++++++ 9 files changed, 312 insertions(+), 1 deletion(-) create mode 100644 src/xps_fitting/diagnostics.py create mode 100644 src/xps_fitting/models.py create mode 100644 src/xps_fitting/optimiser.py create mode 100644 src/xps_fitting/plotting.py create mode 100644 src/xps_fitting/result.py create mode 100644 tests/test_doublet_constraints.py create mode 100644 tests/test_synthetic_fit.py diff --git a/README.md b/README.md index b931ea1..e57083f 100644 --- a/README.md +++ b/README.md @@ -62,3 +62,23 @@ KherveFitting source was not found locally. Its named GL/LA functions can use product or asymmetric conventions, whereas lmfit's `PseudoVoigtModel` is a linear mixture of area-normalised Gaussian and Lorentzian profiles with a shared FWHM. No numerical equivalence is claimed and no Kherve-compatible peak shape is exposed. + +## Fitting and diagnostics + +`fit_spectrum(spectrum, config)` performs staged optimisation: areas, bounded +centres, bounded FWHMs, then optional mixing fractions. Shirley backgrounds are +updated between stages. Seeded multistart fitting detects sensitivity to guesses. +The locally exercised backend is SciPy bounded least squares; `lmfit` remains a +declared project dependency but is not available in the current validation runtime. + +The returned `FitResult` is the Phase 2 contract: ordered energy, raw intensity, +background, labelled component arrays, total fit, residual, fitted parameters, +uncertainties, correlation matrix, fit statistics (RSS, reduced chi-square, AIC, +AICc, BIC, Durbin-Watson, and runs), warnings, configuration, metadata, +convergence details, and software versions. `to_dict()` is JSON-compatible. + +Warnings identify bound hits, high correlations, negligible/unresolved components, +nonphysical backgrounds, convergence problems, and multistart sensitivity. These +diagnostics do not prove chemical correctness. `plot_fit` draws raw points, +background, components, total fit, and an optional residual panel with the binding +energy axis reversed; it can save a diagnostic PNG. diff --git a/src/xps_fitting/__init__.py b/src/xps_fitting/__init__.py index e68bd0e..1713692 100644 --- a/src/xps_fitting/__init__.py +++ b/src/xps_fitting/__init__.py @@ -1,7 +1,8 @@ """Chemistry-constrained XPS fitting tools.""" from .spectrum import Spectrum +from .result import FitResult __version__ = "0.1.0" -__all__ = ["Spectrum", "__version__"] +__all__ = ["FitResult", "Spectrum", "__version__"] diff --git a/src/xps_fitting/diagnostics.py b/src/xps_fitting/diagnostics.py new file mode 100644 index 0000000..5f63c77 --- /dev/null +++ b/src/xps_fitting/diagnostics.py @@ -0,0 +1,19 @@ +"""Residual diagnostics and warning generation.""" + +from __future__ import annotations + +import numpy as np + + +def statistics(residual: np.ndarray, n_parameters: int) -> dict[str, float]: + n = residual.size + rss = float(residual @ residual) + dof = n - n_parameters + safe = max(rss / n, np.finfo(float).tiny) + aic = n * np.log(safe) + 2 * n_parameters + aicc = aic + (2 * n_parameters * (n_parameters + 1) / (n - n_parameters - 1) if n > n_parameters + 1 else np.inf) + bic = n * np.log(safe) + n_parameters * np.log(n) + dw = float(np.diff(residual) @ np.diff(residual) / rss) if rss else 0.0 + signs = residual >= 0 + runs = int(1 + np.count_nonzero(signs[1:] != signs[:-1])) if n else 0 + return {"rss": rss, "reduced_chi_square": rss / dof if dof > 0 else np.nan, "aic": float(aic), "aicc": float(aicc), "bic": float(bic), "durbin_watson": dw, "residual_runs": float(runs), "n_parameters": float(n_parameters), "n_points": float(n)} diff --git a/src/xps_fitting/models.py b/src/xps_fitting/models.py new file mode 100644 index 0000000..31f121e --- /dev/null +++ b/src/xps_fitting/models.py @@ -0,0 +1,41 @@ +"""Numerical evaluation of configured component models.""" + +from __future__ import annotations + +import numpy as np + +from .configuration import PeakConfig +from .lineshapes import gaussian, lorentzian, pseudo_voigt + + +def evaluate_peak(x: np.ndarray, peak: PeakConfig, values: dict[str, float]) -> np.ndarray: + prefix = peak.label + area, centre = values[f"{prefix}.area"], values[f"{prefix}.centre"] + fwhm, fraction = values[f"{prefix}.fwhm"], values[f"{prefix}.fraction"] + if peak.line_shape == "gaussian": + return gaussian(x, area, centre, fwhm) + if peak.line_shape == "lorentzian": + return lorentzian(x, area, centre, fwhm) + if peak.line_shape == "pseudo_voigt": + return pseudo_voigt(x, area, centre, fwhm, fraction) + raise ValueError(f"unsupported fitting line shape: {peak.line_shape}") + + +def resolve_links(peaks: list[PeakConfig], values: dict[str, float]) -> dict[str, float]: + resolved = dict(values) + width_groups: dict[str, float] = {} + fraction_groups: dict[str, float] = {} + for peak in peaks: + if peak.width_group: + key = f"{peak.label}.fwhm" + resolved[key] = width_groups.setdefault(peak.width_group, resolved[key]) + if peak.fraction_group: + key = f"{peak.label}.fraction" + resolved[key] = fraction_groups.setdefault(peak.fraction_group, resolved[key]) + if peak.centre_offset_from: + reference, offset = peak.centre_offset_from + resolved[f"{peak.label}.centre"] = resolved[f"{reference}.centre"] + offset + if peak.area_ratio_to: + reference, ratio = peak.area_ratio_to + resolved[f"{peak.label}.area"] = resolved[f"{reference}.area"] * ratio + return resolved diff --git a/src/xps_fitting/optimiser.py b/src/xps_fitting/optimiser.py new file mode 100644 index 0000000..b5a455f --- /dev/null +++ b/src/xps_fitting/optimiser.py @@ -0,0 +1,118 @@ +"""Deterministic staged constrained fitting.""" + +from __future__ import annotations + +import platform +from dataclasses import replace + +import numpy as np +import scipy +from scipy.optimize import least_squares + +from .backgrounds import linear, shirley +from .configuration import FitConfig, PeakConfig +from .constraints import validate_links +from .diagnostics import statistics +from .models import evaluate_peak, resolve_links +from .result import FitResult +from .spectrum import Spectrum + + +def _parameter_table(peaks: list[PeakConfig]) -> tuple[dict[str, float], dict[str, tuple[float, float]]]: + values, bounds = {}, {} + for peak in peaks: + for name in ("area", "centre", "fwhm", "fraction"): + values[f"{peak.label}.{name}"] = float(getattr(peak, name)) + bounds[f"{peak.label}.{name}"] = tuple(getattr(peak, f"{name}_bounds")) + return values, bounds + + +def _independent_names(peaks: list[PeakConfig], stage: str, release_fraction: bool) -> list[str]: + allowed = {"area"} + if stage in {"centres", "widths", "fractions"}: allowed.add("centre") + if stage in {"widths", "fractions"}: allowed.add("fwhm") + if stage == "fractions" and release_fraction: allowed.add("fraction") + names, widths, fractions = [], set(), set() + for peak in peaks: + for name in ("area", "centre", "fwhm", "fraction"): + if name not in allowed or name in peak.fixed: continue + if name == "area" and peak.area_ratio_to: continue + if name == "centre" and peak.centre_offset_from: continue + if name == "fwhm" and peak.width_group: + if peak.width_group in widths: continue + widths.add(peak.width_group) + if name == "fraction" and peak.fraction_group: + if peak.fraction_group in fractions: continue + fractions.add(peak.fraction_group) + names.append(f"{peak.label}.{name}") + return names + + +def fit_spectrum(spectrum: Spectrum, config: FitConfig, *, backend: str = "scipy") -> FitResult: + """Fit a configured chemical hypothesis using deterministic staged optimisation.""" + if backend != "scipy": + raise ValueError("this installation currently provides the exercised scipy backend") + validate_links(config.peaks) + x, y = spectrum.binding_energy, spectrum.intensity + base_background = linear(x, y[0], y[-1]) if config.background == "linear" else shirley(y) + values, bounds = _parameter_table(config.peaks) + rng = np.random.default_rng(config.random_seed) + best = None + solutions: list[float] = [] + final_names: list[str] = [] + stages = ["areas", "centres", "widths"] + (["fractions"] if config.release_fraction else []) + for start in range(max(1, config.multistart)): + trial = dict(values) + if start: + for peak in config.peaks: + for name in ("area", "centre", "fwhm"): + key = f"{peak.label}.{name}"; lo, hi = bounds[key] + span = (hi - lo) if np.isfinite(hi) else max(abs(trial[key]), 1.0) + trial[key] = float(np.clip(trial[key] + rng.normal(0, 0.08 * span), lo, hi)) + result = None + background = base_background.copy() + for stage in stages: + names = _independent_names(config.peaks, stage, config.release_fraction) + lower = np.array([bounds[name][0] for name in names]); upper = np.array([bounds[name][1] for name in names]) + def residual(vector: np.ndarray) -> np.ndarray: + current = dict(trial); current.update(zip(names, vector)); current = resolve_links(config.peaks, current) + model = sum((evaluate_peak(x, peak, current) for peak in config.peaks), start=np.zeros_like(x)) + return y - background - model + result = least_squares(residual, [trial[name] for name in names], bounds=(lower, upper), loss=config.robust_loss) + trial.update(zip(names, result.x)); trial = resolve_links(config.peaks, trial) + if config.background == "shirley": + components = sum((evaluate_peak(x, peak, trial) for peak in config.peaks), start=np.zeros_like(x)) + background = shirley(y - components) + components * 0 + assert result is not None + score = float(result.fun @ result.fun); solutions.append(score) + if best is None or score < best[0]: best = (score, result, trial, background, names) + assert best is not None + _, raw_result, fitted, background, final_names = best + components = {peak.label: evaluate_peak(x, peak, fitted) for peak in config.peaks} + total = background + sum(components.values(), start=np.zeros_like(x)); residual = y - total + stats = statistics(residual, len(final_names)) + uncertainties: dict[str, float | None] = {name: None for name in fitted}; correlations: dict[str, dict[str, float]] = {} + warnings: list[str] = [] + if raw_result.jac.size and raw_result.jac.shape[0] > raw_result.jac.shape[1]: + try: + covariance = np.linalg.inv(raw_result.jac.T @ raw_result.jac) * stats["reduced_chi_square"] + stderr = np.sqrt(np.diag(covariance)); corr = covariance / np.outer(stderr, stderr) + uncertainties.update(dict(zip(final_names, map(float, stderr)))) + correlations = {a: {b: float(corr[i, j]) for j, b in enumerate(final_names)} for i, a in enumerate(final_names)} + if np.any(np.abs(corr - np.eye(len(corr))) > 0.95): warnings.append("Parameter correlation exceeds 0.95.") + except np.linalg.LinAlgError: warnings.append("Covariance matrix is singular; uncertainties are unavailable.") + for peak in config.peaks: + for name in ("area", "centre", "fwhm", "fraction"): + key = f"{peak.label}.{name}"; lo, hi = bounds[key]; value = fitted[key] + scale = max(1.0, abs(value), abs(hi - lo) if np.isfinite(hi) else 1.0) + if abs(value - lo) <= 1e-5 * scale or (np.isfinite(hi) and abs(hi - value) <= 1e-5 * scale): warnings.append(f"{key} is at or near a bound.") + if fitted[f"{peak.label}.area"] < 0.005 * sum(p.area for p in config.peaks): warnings.append(f"{peak.label} is a negligible component.") + for index, first in enumerate(config.peaks): + for second in config.peaks[index + 1:]: + separation = abs(fitted[f"{first.label}.centre"] - fitted[f"{second.label}.centre"]) + mean_width = (fitted[f"{first.label}.fwhm"] + fitted[f"{second.label}.fwhm"]) / 2 + if separation < 0.5 * mean_width: warnings.append(f"{first.label} and {second.label} are unresolved by the 0.5 FWHM criterion.") + if not raw_result.success: warnings.append(f"Convergence warning: {raw_result.message}") + if np.any(background < 0): warnings.append("Background contains negative values.") + if len(solutions) > 1 and (max(solutions) - min(solutions)) / max(min(solutions), 1e-30) > 0.05: warnings.append("Multiple initialisations produced materially different solutions.") + return FitResult(x, y, background, components, total, residual, fitted, uncertainties, correlations, stats, warnings, config.to_dict(), dict(spectrum.metadata), {"success": bool(raw_result.success), "message": raw_result.message, "evaluations": raw_result.nfev, "backend": backend, "multistart_rss": solutions}, {"xps_fitting": "0.1.0", "python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__}) diff --git a/src/xps_fitting/plotting.py b/src/xps_fitting/plotting.py new file mode 100644 index 0000000..5994fd4 --- /dev/null +++ b/src/xps_fitting/plotting.py @@ -0,0 +1,26 @@ +"""Clean diagnostic plotting (publication themes are Phase 2).""" + +from __future__ import annotations + +from pathlib import Path + +import matplotlib.pyplot as plt + +from .result import FitResult + + +def plot_fit(result: FitResult, path: str | Path | None = None, *, residual_panel: bool = True): + figure, axes = plt.subplots(2 if residual_panel else 1, 1, sharex=True, gridspec_kw={"height_ratios": [3, 1]} if residual_panel else None) + main = axes[0] if residual_panel else axes + main.plot(result.energy, result.raw_intensity, ".", color="0.35", label="data") + main.plot(result.energy, result.background, "--", color="0.5", label="background") + for label, curve in result.components.items(): main.plot(result.energy, result.background + curve, lw=1, label=label) + main.plot(result.energy, result.total_fit, color="black", label="total fit") + main.legend(); main.set_ylabel("Intensity"); main.invert_xaxis() + if residual_panel: + axes[1].axhline(0, color="0.5", lw=0.8); axes[1].plot(result.energy, result.residual, color="black") + axes[1].set_ylabel("Residual"); axes[1].set_xlabel("Binding energy (eV)") + else: main.set_xlabel("Binding energy (eV)") + figure.tight_layout() + if path is not None: figure.savefig(path, dpi=200) + return figure diff --git a/src/xps_fitting/result.py b/src/xps_fitting/result.py new file mode 100644 index 0000000..f060302 --- /dev/null +++ b/src/xps_fitting/result.py @@ -0,0 +1,42 @@ +"""Stable numerical result contract for fitting and Phase 2 plotting.""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass, field +from typing import Any + +import numpy as np + + +def _jsonable(value: Any) -> Any: + if isinstance(value, np.ndarray): + return value.tolist() + if isinstance(value, dict): + return {str(key): _jsonable(item) for key, item in value.items()} + if isinstance(value, (list, tuple)): + return [_jsonable(item) for item in value] + if isinstance(value, np.generic): + return value.item() + return value + + +@dataclass +class FitResult: + energy: np.ndarray + raw_intensity: np.ndarray + background: np.ndarray + components: dict[str, np.ndarray] + total_fit: np.ndarray + residual: np.ndarray + fitted_parameters: dict[str, float] + parameter_uncertainties: dict[str, float | None] = field(default_factory=dict) + correlation_matrix: dict[str, dict[str, float]] = field(default_factory=dict) + fit_statistics: dict[str, float] = field(default_factory=dict) + warnings: list[str] = field(default_factory=list) + configuration: dict[str, Any] = field(default_factory=dict) + metadata: dict[str, Any] = field(default_factory=dict) + convergence: dict[str, Any] = field(default_factory=dict) + software_versions: dict[str, str] = field(default_factory=dict) + + def to_dict(self) -> dict[str, Any]: + return _jsonable(asdict(self)) diff --git a/tests/test_doublet_constraints.py b/tests/test_doublet_constraints.py new file mode 100644 index 0000000..4b1fb54 --- /dev/null +++ b/tests/test_doublet_constraints.py @@ -0,0 +1,20 @@ +import numpy as np + +from xps_fitting.configuration import FitConfig +from xps_fitting.constraints import cl2p_doublet +from xps_fitting.lineshapes import pseudo_voigt +from xps_fitting.optimiser import fit_spectrum +from xps_fitting.spectrum import Spectrum + + +def test_fitted_doublet_links_are_exact() -> None: + x = np.linspace(195, 207, 401) + y = 10 + pseudo_voigt(x, 900, 200, 1.2, 0.4) + pseudo_voigt(x, 450, 201.6, 1.2, 0.4) + result = fit_spectrum(Spectrum(x, y), FitConfig("cl", "Cl 2p", cl2p_doublet("Cl", 199.8, 800, fraction=0.4))) + p = result.fitted_parameters + assert p["Cl_2p1/2.centre"] == pytest.approx(p["Cl_2p3/2.centre"] + 1.6) + assert p["Cl_2p1/2.area"] == pytest.approx(p["Cl_2p3/2.area"] / 2) + assert p["Cl_2p1/2.fwhm"] == p["Cl_2p3/2.fwhm"] + + +import pytest diff --git a/tests/test_synthetic_fit.py b/tests/test_synthetic_fit.py new file mode 100644 index 0000000..338d6d7 --- /dev/null +++ b/tests/test_synthetic_fit.py @@ -0,0 +1,24 @@ +import numpy as np + +from xps_fitting.configuration import FitConfig, PeakConfig +from xps_fitting.lineshapes import pseudo_voigt +from xps_fitting.optimiser import fit_spectrum +from xps_fitting.plotting import plot_fit +from xps_fitting.spectrum import Spectrum + + +def test_synthetic_recovery_is_deterministic(tmp_path) -> None: + x = np.linspace(280, 292, 481) + rng = np.random.default_rng(5) + y = 20 + pseudo_voigt(x, 1000, 285.0, 1.3, 0.3) + rng.normal(0, 0.15, x.size) + config = FitConfig("one", "C 1s", [PeakConfig("C-C", 284.8, (284.5, 285.5), 900, fwhm=1.5, fwhm_bounds=(1, 2), fraction=0.3, fixed=("fraction",))], multistart=2, random_seed=9) + first = fit_spectrum(Spectrum(x, y, region="C 1s"), config) + second = fit_spectrum(Spectrum(x, y, region="C 1s"), config) + assert first.fitted_parameters["C-C.centre"] == pytest.approx(285, abs=0.01) + assert first.fitted_parameters["C-C.area"] == pytest.approx(1000, rel=0.02) + assert first.fit_statistics["rss"] == second.fit_statistics["rss"] + output = tmp_path / "diagnostic.png"; plot_fit(first, output) + assert output.stat().st_size > 0 and first.to_dict()["energy"][0] == 280 + + +import pytest From d43752c16697ed62b5bd50e21a925dcefa76fa5d Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:56:29 +0800 Subject: [PATCH 05/75] feat: add chemistry-informed models and model comparison --- README.md | 20 ++++++++++++++++++++ configs/pdi_h_cooh_c1s_4.json | 9 +++++++++ configs/pdi_h_cooh_c1s_5.json | 10 ++++++++++ examples/fit_cl2p_doublet.py | 19 +++++++++++++++++++ examples/fit_pdi_h_cooh_c1s.py | 20 ++++++++++++++++++++ src/xps_fitting/configuration.py | 10 ++++++++++ src/xps_fitting/model_comparison.py | 27 +++++++++++++++++++++++++++ src/xps_fitting/optimiser.py | 11 +++++++++-- tests/test_model_comparison.py | 18 ++++++++++++++++++ 9 files changed, 142 insertions(+), 2 deletions(-) create mode 100644 configs/pdi_h_cooh_c1s_4.json create mode 100644 configs/pdi_h_cooh_c1s_5.json create mode 100644 examples/fit_cl2p_doublet.py create mode 100644 examples/fit_pdi_h_cooh_c1s.py create mode 100644 src/xps_fitting/model_comparison.py create mode 100644 tests/test_model_comparison.py diff --git a/README.md b/README.md index e57083f..546decb 100644 --- a/README.md +++ b/README.md @@ -82,3 +82,23 @@ nonphysical backgrounds, convergence problems, and multistart sensitivity. These diagnostics do not prove chemical correctness. `plot_fit` draws raw points, background, components, total fit, and an optional residual panel with the binding energy axis reversed; it can save a diagnostic PNG. + +## Candidate models + +JSON configuration files in `configs/` define the four- and five-component +PDI-H-COOH C 1s hypotheses. The four assignments are aromatic C-C/C=C (284.65 eV), +C-N/C-Cl (285.85 eV), imide N-C=O (287.90 eV), and acid O-C=O (289.15 eV). The +five-component alternative adds a broader pi-pi* satellite near 290.70 eV. There +is deliberately no default 284.10 eV component, and ordinary widths have separate +1.1–2.0 eV bounds rather than an exact global equality. Config objects also support +shared-width groups, a width soft-penalty setting, and relative satellite offsets. + +`compare_models` returns every `FitResult`; `comparison_table` reports AIC/AICc/BIC, +RSS, residual diagnostics, warnings, and multistart stability. Resolution and +meaningfulness warnings must be reviewed alongside chemical knowledge: a lower +information criterion is not a declaration of chemical truth. + +```bash +MPLBACKEND=Agg python examples/fit_pdi_h_cooh_c1s.py +MPLBACKEND=Agg python examples/fit_cl2p_doublet.py +``` diff --git a/configs/pdi_h_cooh_c1s_4.json b/configs/pdi_h_cooh_c1s_4.json new file mode 100644 index 0000000..e51026f --- /dev/null +++ b/configs/pdi_h_cooh_c1s_4.json @@ -0,0 +1,9 @@ +{ + "name": "C1s_4", "region": "C 1s", "background": "linear", "multistart": 2, "random_seed": 42, + "peaks": [ + {"label": "aromatic_C-C_C=C", "centre": 284.65, "centre_bounds": [284.35, 284.95], "area": 1000, "fwhm": 1.4, "fwhm_bounds": [1.1, 2.0]}, + {"label": "C-N_C-Cl", "centre": 285.85, "centre_bounds": [285.3, 286.4], "area": 500, "fwhm": 1.4, "fwhm_bounds": [1.1, 2.0]}, + {"label": "imide_N-C=O", "centre": 287.90, "centre_bounds": [287.4, 288.4], "area": 400, "fwhm": 1.5, "fwhm_bounds": [1.1, 2.0]}, + {"label": "acid_O-C=O", "centre": 289.15, "centre_bounds": [288.7, 289.7], "area": 250, "fwhm": 1.5, "fwhm_bounds": [1.1, 2.0]} + ] +} diff --git a/configs/pdi_h_cooh_c1s_5.json b/configs/pdi_h_cooh_c1s_5.json new file mode 100644 index 0000000..1d7c773 --- /dev/null +++ b/configs/pdi_h_cooh_c1s_5.json @@ -0,0 +1,10 @@ +{ + "name": "C1s_5", "region": "C 1s", "background": "linear", "multistart": 2, "random_seed": 42, + "peaks": [ + {"label": "aromatic_C-C_C=C", "centre": 284.65, "centre_bounds": [284.35, 284.95], "area": 1000, "fwhm": 1.4, "fwhm_bounds": [1.1, 2.0]}, + {"label": "C-N_C-Cl", "centre": 285.85, "centre_bounds": [285.3, 286.4], "area": 500, "fwhm": 1.4, "fwhm_bounds": [1.1, 2.0]}, + {"label": "imide_N-C=O", "centre": 287.90, "centre_bounds": [287.4, 288.4], "area": 400, "fwhm": 1.5, "fwhm_bounds": [1.1, 2.0]}, + {"label": "acid_O-C=O", "centre": 289.15, "centre_bounds": [288.7, 289.7], "area": 250, "fwhm": 1.5, "fwhm_bounds": [1.1, 2.0]}, + {"label": "pi-pi_star", "centre": 290.70, "centre_bounds": [290.1, 291.8], "area": 100, "fwhm": 2.2, "fwhm_bounds": [1.5, 3.5], "centre_offset_from": ["aromatic_C-C_C=C", 6.05]} + ] +} diff --git a/examples/fit_cl2p_doublet.py b/examples/fit_cl2p_doublet.py new file mode 100644 index 0000000..67c4e78 --- /dev/null +++ b/examples/fit_cl2p_doublet.py @@ -0,0 +1,19 @@ +"""Fit a constrained synthetic Cl 2p doublet.""" + +from pathlib import Path + +import numpy as np + +from xps_fitting.configuration import FitConfig +from xps_fitting.constraints import cl2p_doublet +from xps_fitting.lineshapes import pseudo_voigt +from xps_fitting.optimiser import fit_spectrum +from xps_fitting.plotting import plot_fit +from xps_fitting.spectrum import Spectrum + +x = np.linspace(195, 207, 481) +y = 10 + pseudo_voigt(x, 1000, 200, 1.2, 0.5) + pseudo_voigt(x, 500, 201.6, 1.2, 0.5) +result = fit_spectrum(Spectrum(x, y, region="Cl 2p"), FitConfig("Cl_doublet", "Cl 2p", cl2p_doublet("Cl", 199.8, 900))) +print(result.fit_statistics) +output = Path(__file__).resolve().parents[1] / "outputs"; output.mkdir(exist_ok=True) +plot_fit(result, output / "cl2p_diagnostic.png") diff --git a/examples/fit_pdi_h_cooh_c1s.py b/examples/fit_pdi_h_cooh_c1s.py new file mode 100644 index 0000000..8c77a42 --- /dev/null +++ b/examples/fit_pdi_h_cooh_c1s.py @@ -0,0 +1,20 @@ +"""Fit and compare synthetic PDI-H-COOH C 1s candidate hypotheses.""" + +from pathlib import Path + +import numpy as np + +from xps_fitting.configuration import load_config +from xps_fitting.lineshapes import pseudo_voigt +from xps_fitting.model_comparison import compare_models, comparison_table +from xps_fitting.plotting import plot_fit +from xps_fitting.spectrum import Spectrum + +root = Path(__file__).resolve().parents[1] +x = np.linspace(280, 294, 561) +y = 25 + sum(pseudo_voigt(x, area, centre, width, 0.5) for area, centre, width in [(1200, 284.65, 1.35), (500, 285.85, 1.45), (400, 287.9, 1.5), (250, 289.15, 1.55), (120, 290.7, 2.2)]) +configs = [load_config(root / "configs" / f"pdi_h_cooh_c1s_{count}.json") for count in (4, 5)] +results = compare_models(Spectrum(x, y, region="C 1s", sample_name="synthetic PDI-H-COOH"), configs) +print(comparison_table(results)) +output = root / "outputs"; output.mkdir(exist_ok=True) +plot_fit(results["C1s_5"], output / "pdi_h_cooh_c1s_diagnostic.png") diff --git a/src/xps_fitting/configuration.py b/src/xps_fitting/configuration.py index c746fc1..71dd3f4 100644 --- a/src/xps_fitting/configuration.py +++ b/src/xps_fitting/configuration.py @@ -3,6 +3,8 @@ from __future__ import annotations from dataclasses import asdict, dataclass, field +import json +from pathlib import Path from typing import Any @@ -48,3 +50,11 @@ class FitConfig: def to_dict(self) -> dict[str, Any]: return asdict(self) + + +def load_config(path: str | Path) -> FitConfig: + """Load a JSON candidate-model configuration.""" + with Path(path).open(encoding="utf-8") as stream: + data = json.load(stream) + data["peaks"] = [PeakConfig(**peak) for peak in data["peaks"]] + return FitConfig(**data) diff --git a/src/xps_fitting/model_comparison.py b/src/xps_fitting/model_comparison.py new file mode 100644 index 0000000..9ebc2a3 --- /dev/null +++ b/src/xps_fitting/model_comparison.py @@ -0,0 +1,27 @@ +"""Compare candidate hypotheses without asserting chemical correctness.""" + +from __future__ import annotations + +from .configuration import FitConfig +from .optimiser import fit_spectrum +from .result import FitResult +from .spectrum import Spectrum + + +def compare_models(spectrum: Spectrum, configurations: list[FitConfig]) -> dict[str, FitResult]: + if len({config.name for config in configurations}) != len(configurations): + raise ValueError("candidate model names must be unique") + return {config.name: fit_spectrum(spectrum, config) for config in configurations} + + +def comparison_table(results: dict[str, FitResult]) -> list[dict[str, object]]: + table = [] + for name, result in results.items(): + table.append({ + "model": name, + **{key: result.fit_statistics[key] for key in ("aic", "aicc", "bic", "rss", "durbin_watson", "residual_runs")}, + "warning_count": len(result.warnings), + "warnings": list(result.warnings), + "multistart_rss": result.convergence.get("multistart_rss", []), + }) + return table diff --git a/src/xps_fitting/optimiser.py b/src/xps_fitting/optimiser.py index b5a455f..4140586 100644 --- a/src/xps_fitting/optimiser.py +++ b/src/xps_fitting/optimiser.py @@ -77,14 +77,21 @@ def fit_spectrum(spectrum: Spectrum, config: FitConfig, *, backend: str = "scipy def residual(vector: np.ndarray) -> np.ndarray: current = dict(trial); current.update(zip(names, vector)); current = resolve_links(config.peaks, current) model = sum((evaluate_peak(x, peak, current) for peak in config.peaks), start=np.zeros_like(x)) - return y - background - model + data_residual = y - background - model + if config.width_penalty > 0 and len(config.peaks) > 1: + widths = np.array([current[f"{peak.label}.fwhm"] for peak in config.peaks]) + penalty = np.sqrt(config.width_penalty) * (widths[1:] - widths[:-1]) + return np.concatenate((data_residual, penalty)) + return data_residual result = least_squares(residual, [trial[name] for name in names], bounds=(lower, upper), loss=config.robust_loss) trial.update(zip(names, result.x)); trial = resolve_links(config.peaks, trial) if config.background == "shirley": components = sum((evaluate_peak(x, peak, trial) for peak in config.peaks), start=np.zeros_like(x)) background = shirley(y - components) + components * 0 assert result is not None - score = float(result.fun @ result.fun); solutions.append(score) + current = resolve_links(config.peaks, trial) + physical_residual = y - background - sum((evaluate_peak(x, peak, current) for peak in config.peaks), start=np.zeros_like(x)) + score = float(physical_residual @ physical_residual); solutions.append(score) if best is None or score < best[0]: best = (score, result, trial, background, names) assert best is not None _, raw_result, fitted, background, final_names = best diff --git a/tests/test_model_comparison.py b/tests/test_model_comparison.py new file mode 100644 index 0000000..1c276c4 --- /dev/null +++ b/tests/test_model_comparison.py @@ -0,0 +1,18 @@ +from pathlib import Path + +import numpy as np + +from xps_fitting.configuration import load_config +from xps_fitting.lineshapes import pseudo_voigt +from xps_fitting.model_comparison import compare_models, comparison_table +from xps_fitting.spectrum import Spectrum + + +def test_five_component_signal_supports_five_component_candidate() -> None: + root = Path(__file__).parents[1] + x = np.linspace(280, 294, 401) + y = 20 + sum(pseudo_voigt(x, a, c, w, 0.5) for a, c, w in [(1000, 284.65, 1.4), (500, 285.85, 1.4), (400, 287.9, 1.5), (250, 289.15, 1.5), (150, 290.7, 2.2)]) + configs = [load_config(root / "configs" / f"pdi_h_cooh_c1s_{n}.json") for n in (4, 5)] + table = comparison_table(compare_models(Spectrum(x, y), configs)) + assert {row["model"] for row in table} == {"C1s_4", "C1s_5"} + assert next(row for row in table if row["model"] == "C1s_5")["bic"] < next(row for row in table if row["model"] == "C1s_4")["bic"] From 24a478e5457db54a8b5fc6739fbada63d0d2c5ed Mon Sep 17 00:00:00 2001 From: Squiddy Date: Mon, 20 Jul 2026 23:58:53 +0800 Subject: [PATCH 06/75] feat: complete phase one fitting and export workflow --- CHANGELOG.md | 9 +++++++ README.md | 33 +++++++++++++++++++++++++ docs/fitresult.schema.json | 18 ++++++++++++++ docs/fitresult_contract.md | 11 +++++++++ docs/fitting_methodology.md | 15 ++++++++++++ src/xps_fitting/export.py | 45 +++++++++++++++++++++++++++++++++++ src/xps_fitting/global_fit.py | 40 +++++++++++++++++++++++++++++++ src/xps_fitting/optimiser.py | 5 +++- src/xps_fitting/result.py | 5 +++- tests/test_export.py | 20 ++++++++++++++++ tests/test_global_fit.py | 15 ++++++++++++ 11 files changed, 214 insertions(+), 2 deletions(-) create mode 100644 CHANGELOG.md create mode 100644 docs/fitresult.schema.json create mode 100644 docs/fitresult_contract.md create mode 100644 docs/fitting_methodology.md create mode 100644 src/xps_fitting/export.py create mode 100644 src/xps_fitting/global_fit.py create mode 100644 tests/test_export.py create mode 100644 tests/test_global_fit.py diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..83a74e0 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,9 @@ +# Changelog + +## Unreleased — Phase 1 + +- Added validated CSV, XLSX, DataFrame, and optional VGD spectrum input. +- Added area-normalised line shapes, backgrounds, chemistry constraints, staged + deterministic fitting, diagnostics, model comparison, and diagnostic plotting. +- Added PDI-H-COOH C 1s and Cl 2p examples, export bundles, a stable `FitResult` + contract, and limited two-pass shared-shape fitting groundwork. diff --git a/README.md b/README.md index 546decb..f99748b 100644 --- a/README.md +++ b/README.md @@ -102,3 +102,36 @@ information criterion is not a declaration of chemical truth. MPLBACKEND=Agg python examples/fit_pdi_h_cooh_c1s.py MPLBACKEND=Agg python examples/fit_cl2p_doublet.py ``` + +## Exports and linked fitting + +`export_result(result, directory)` writes a curve CSV, an XLSX workbook with curves, +parameters, statistics, warnings and metadata, a JSON parameter/configuration +summary, a Markdown report, and a diagnostic PNG. Generated files belong in the +ignored `outputs/` directory. The complete result API and schema are documented in +[`docs/fitresult_contract.md`](docs/fitresult_contract.md). + +`fit_shared_shapes` provides limited global-fitting groundwork: it performs separate +first-pass fits, averages selected FWHM/fraction values, and fixes those consensus +values during a second pass. Areas, centres, and intensity remain sample-specific; +configured relative offsets stay exact. This is not simultaneous optimisation, +does not estimate joint covariance, and does not yet fit a sample-wide charging +shift or bounded per-sample deviations. + +## Reproducibility and development + +Use explicit candidate configurations, preserve raw intensity, record crop and +acquisition metadata, set the multistart seed, inspect warnings/residuals, and export +the result bundle. Run the full validation and examples with: + +```bash +MPLBACKEND=Agg pytest +PYTHONPATH=src MPLBACKEND=Agg python examples/load_spectrum.py +PYTHONPATH=src MPLBACKEND=Agg python examples/fit_pdi_h_cooh_c1s.py +PYTHONPATH=src MPLBACKEND=Agg python examples/fit_cl2p_doublet.py +``` + +Current limitations include an unexercised lmfit dependency in this environment, +no Tougaard/asymmetric lines, no uncertainty weighting, and approximate rather than +joint global fitting. See [the methodology](docs/fitting_methodology.md). Phase 2 +will build publication-quality styling strictly on the `FitResult` arrays. diff --git a/docs/fitresult.schema.json b/docs/fitresult.schema.json new file mode 100644 index 0000000..630b5bd --- /dev/null +++ b/docs/fitresult.schema.json @@ -0,0 +1,18 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "title": "xps_fitting.FitResult", + "type": "object", + "required": ["energy", "raw_intensity", "background", "components", "total_fit", "residual", "fitted_parameters", "parameter_uncertainties", "correlation_matrix", "fit_statistics", "warnings", "configuration", "metadata", "convergence", "software_versions"], + "properties": { + "energy": {"type": "array", "items": {"type": "number"}}, + "raw_intensity": {"type": "array", "items": {"type": "number"}}, + "background": {"type": "array", "items": {"type": "number"}}, + "components": {"type": "object", "additionalProperties": {"type": "array", "items": {"type": "number"}}}, + "total_fit": {"type": "array", "items": {"type": "number"}}, + "residual": {"type": "array", "items": {"type": "number"}}, + "fitted_parameters": {"type": "object", "additionalProperties": {"type": "number"}}, + "parameter_uncertainties": {"type": "object"}, "correlation_matrix": {"type": "object"}, + "fit_statistics": {"type": "object"}, "warnings": {"type": "array", "items": {"type": "string"}}, + "configuration": {"type": "object"}, "metadata": {"type": "object"}, "convergence": {"type": "object"}, "software_versions": {"type": "object"} + } +} diff --git a/docs/fitresult_contract.md b/docs/fitresult_contract.md new file mode 100644 index 0000000..a030ef9 --- /dev/null +++ b/docs/fitresult_contract.md @@ -0,0 +1,11 @@ +# FitResult contract + +`xps_fitting.result.FitResult` is the Phase 2 numerical interface. Its curve arrays +share one ascending energy grid. `components` contains background-free labelled peak +arrays; `total_fit` equals background plus their sum and `residual` equals raw minus +total. Parameters use `