From 73a5e3efc74567140409c135a75fd4f6a17fd80a Mon Sep 17 00:00:00 2001 From: vedina Date: Mon, 5 Oct 2026 19:58:39 +0300 Subject: [PATCH 1/2] Add nexus_calibration: export a ramanchada2 calibration as NXraman NeXus Moved from ramanchada2.protocols.calibration.serialization so ramanchada2 does not depend on pyambit or nexusformat. The module takes a derived CalibrationModel (and optional YCalibrationComponent / calibrant spectra), writes the spectra through spe2ambit / NXRamanProtocolApplication and adds NXcalibration groups (curve, fit parameters, matched-peak anchors, JSON model) via nexusformat. ramanchada2 is imported lazily; the tests skip when it is not installed. Co-Authored-By: Claude Sonnet 5.5 --- src/pyambit/nexus_calibration.py | 443 ++++++++++++++++++ .../nexus_models/nexus_calibration_test.py | 296 ++++++++++++ 2 files changed, 739 insertions(+) create mode 100644 src/pyambit/nexus_calibration.py create mode 100644 tests/pyambit/nexus_models/nexus_calibration_test.py diff --git a/src/pyambit/nexus_calibration.py b/src/pyambit/nexus_calibration.py new file mode 100644 index 0000000..414d115 --- /dev/null +++ b/src/pyambit/nexus_calibration.py @@ -0,0 +1,443 @@ +"""Export a ramanchada2 calibration workflow as a self-describing NeXus/HDF5 file. + +Moved here from ``ramanchada2.protocols.calibration.serialization`` so that ramanchada2 +does not depend on pyambit or nexusformat. This module needs ramanchada2 (a derived +``CalibrationModel`` / ``YCalibrationComponent`` and ``Spectrum`` objects), which is +imported lazily; ramanchada2 never imports pyambit, so there is no dependency cycle. +""" +import json + +import numpy as np + +SI_REF_CM1 = 520.45 + + +def _calibration_curve(calmodel, spectral_range, npoints): + """Sample the model over ``spectral_range`` (cm-1): uncalibrated -> calibrated.""" + from ramanchada2.spectrum import Spectrum + + lo, hi = float(min(spectral_range)), float(max(spectral_range)) + grid = np.linspace(lo, hi, int(npoints)) + spe = Spectrum(x=grid, y=np.ones_like(grid)) + calibrated = calmodel.apply_calibration_x(spe, spe_units="cm-1").x + return grid, np.asarray(calibrated, dtype=float) + + +def _laser_zero_info(calmodel): + """Si peak position (nm and implied cm-1 on the nominal axis) and calibrated laser wl.""" + from ramanchada2.protocols.calibration.xcalibration import LazerZeroingComponent + + for c in calmodel.components: + if isinstance(c, LazerZeroingComponent): + zero_nm = float(c.model) + si_ref = float(list(c.ref.keys())[0]) if c.ref else SI_REF_CM1 + laser_nm = 1e7 / (1e7 / zero_nm + si_ref) + return { + "si_peak_nm": zero_nm, + "si_reference_cm1": si_ref, + "calibrated_laser_wl_nm": laser_nm, + } + return {} + + +def _h5_clean(value): + """Coerce a to_dict()-sourced value into something h5py's create_dataset accepts. + + Plain numpy string arrays (numpy's default fixed-width unicode dtype, e.g. ' 1 else "calibration_x" + inst[name] = _nx_calibration_group(comp_dict, curve=curve) + + if ycal_component is not None: + y_dict = ycal_component.to_dict() + inst["calibration_y"] = _nx_calibration_group(y_dict) + # Plottable y-calibration curve, the analogue of calibration_curve for x -- the + # gap that made the earlier h5py-only version look like intensity calibration was + # silently dropped (it wasn't; it just had no plottable representation). + if "calibration_curve_y_measured" in (spe_names := [s[3] for s in spe_list]): + y_curve = nx.NXdata() + idx = spe_names.index("calibration_curve_y_measured") + y_curve["calibrated_cm1"] = np.asarray(spe_list[idx][0], dtype=float) + y_curve["intensity_factor"] = np.asarray(spe_list[idx][1], dtype=float) + y_curve.attrs["signal"] = "intensity_factor" + y_curve.attrs["axes"] = ["calibrated_cm1"] + y_curve.attrs["description"] = ("Measured SRM response, resampled onto the " + "certificate's declared range.") + entry["calibration_curve_y"] = y_curve + + nx_root.save(filename, mode="w") + return filename + + +def _nx_calibration_group(component_dict, curve=None): + """Build one calibration component as an NXcalibration nexusformat.nexus.tree group, + for composing into the shared NXroot alongside the pyambit-written entry (the whole + file is saved in one nx_root.save() call — see export_nexus_calibration). ``curve`` + is an optional (original_axis, calibrated_axis) array pair — real datasets, not the + scalar-typed fields pyambit's generated NXCalibration pydantic model exposes (that + codegen has no NX_FLOAT[rank] support yet; see docs/nexus_export_plan.md). + calibration_parameters is a real NXparameters container with its actual numeric + children, not the near-empty base-class stub; calibration_object/NXnote is kept + deliberately small — only fields with no typed home yet (see + _residual_component_fields) — so the typed datasets are the single source of truth + for the reconstructable numbers rather than a duplicate copy. + """ + import nexusformat.nexus.tree as nx + + cal = nx.NXgroup() + cal.nxclass = "NXcalibration" + cal.attrs["description"] = "ramanchada2 " + component_dict.get("type", "calibration") + cal.attrs["physical_quantity"] = ( + "relative intensity" if component_dict.get("type") == "YCalibrationComponent" + else "wavenumber") + cal.attrs["applied"] = bool(component_dict.get("enabled", True)) + cal.attrs["fit_formula_description"] = _fit_formula_description(component_dict) + + if curve is not None: + original_axis, calibrated_axis = curve + cal["original_axis"] = np.asarray(original_axis, dtype=float) + cal["calibrated_axis"] = np.asarray(calibrated_axis, dtype=float) + + params = _component_parameters(component_dict) + if params: + pgroup = nx.NXgroup() + pgroup.nxclass = "NXparameters" + for key, value in params.items(): + if value is None: + continue + pgroup[key] = _h5_clean(value) + cal["calibration_parameters"] = pgroup + + note = nx.NXnote() + note["type"] = "application/json" + note["data"] = json.dumps(_residual_component_fields(component_dict)) + cal["calibration_object"] = note + + anchors = component_dict.get("anchors") + if anchors: + agroup = nx.NXdata() + agroup.attrs["description"] = "Matched calibrant peaks used to derive the model." + arr = np.asarray(anchors, dtype=float) + agroup["measured"] = arr[:, 0] + agroup["reference"] = arr[:, 1] + agroup["inlier"] = arr[:, 2].astype(bool) + cal["anchors"] = agroup + + return cal diff --git a/tests/pyambit/nexus_models/nexus_calibration_test.py b/tests/pyambit/nexus_models/nexus_calibration_test.py new file mode 100644 index 0000000..b030988 --- /dev/null +++ b/tests/pyambit/nexus_models/nexus_calibration_test.py @@ -0,0 +1,296 @@ +#!/usr/bin/env python3 +"""Tests for export_nexus_calibration: the calibration workflow as a NeXus/HDF5 file. + +Reuses the same real-data fixtures as test_calmodel_serialization.py (from_test_spe / +calibration_model_factory), so these exercise a genuine derived model, not a synthetic +stand-in. +""" + +import json + +import h5py +import numpy as np +import pytest + +pytest.importorskip("ramanchada2") + +# the shared fixtures use CalibrationModel.calibration_model_factory (deprecated in rc2) +pytestmark = pytest.mark.filterwarnings("ignore::DeprecationWarning") + +import ramanchada2.misc.constants as rc2const # noqa: E402 +from ramanchada2.protocols.calibration.calibration_model import CalibrationModel # noqa: E402 +from pyambit.nexus_calibration import export_nexus_calibration # noqa: E402 +from ramanchada2.protocols.calibration.ycalibration import ( # noqa: E402 + YCalibrationCertificate, + YCalibrationComponent, +) +from ramanchada2.spectrum import from_test_spe # noqa: E402 + +LASER_WL = 785 +_SPE_KW = dict(provider=["ICV"], device=["BWtek"], OP=["100"], laser_wl=[str(LASER_WL)]) + + +def _load(sample): + return from_test_spe(sample=[sample], **_SPE_KW) + + +@pytest.fixture(scope="module") +def spe_neon(): + return _load("Neon").trim_axes(method="x-axis", boundaries=(100, 3500)) + + +@pytest.fixture(scope="module") +def spe_sil(): + return _load("S0B").trim_axes(method="x-axis", boundaries=(520.45 - 50, 520.45 + 50)) + + +@pytest.fixture(scope="module", params=["poly", "pchipinverse"]) +def calmodel(request, spe_neon, spe_sil): + spe_neon_bl = spe_neon.subtract_baseline_rc1_snip(niter=40) + spe_sil_bl = spe_sil.subtract_baseline_rc1_snip(niter=40) + return CalibrationModel.calibration_model_factory( + LASER_WL, + spe_neon_bl, + spe_sil_bl, + neon_wl=rc2const.NEON_WL[LASER_WL], + find_kw={"wlen": 200, "width": 1}, + fit_peaks_kw={}, + should_fit=False, + interpolator_method=request.param, + ) + + +@pytest.fixture(scope="module") +def ycal(): + cert = YCalibrationCertificate.load(wavelength=LASER_WL, key="NIST785_SRM2241") + spe_srm = _load("NIST785_SRM2241") + return YCalibrationComponent(LASER_WL, reference_spe_xcalibrated=spe_srm, certificate=cert) + + +@pytest.fixture(scope="module") +def spe_pst(): + return _load("PST").trim_axes(method="x-axis", boundaries=(100, 3500)) + + +def _write(calmodel, tmp_path, **kw): + path = str(tmp_path / "calibration.nxs") + export_nexus_calibration(calmodel, path, spectral_range=(150.0, 3400.0), npoints=101, **kw) + return path + + +def test_build_measurement_papp_entry_is_discoverable(tmp_path): + """Regression: pyambit's to_nexus() writes the entry under a leading-slash path + derived from (provider, papp.uuid); it's reachable via nx_root.entries (the pattern + pyambit's own test suite uses), not .keys(). Also: spe2ambit's per-effect group name + comes from its `sample` argument (routed to spe2effect's nx_name) -- passing the same + generic label for every spectrum makes every effect group indistinguishable except by + position, so _build_measurement_papp must pass each spectrum's own nx_name there.""" + from pyambit.nexus_calibration import _build_measurement_papp + import nexusformat.nexus.tree as nx + + spe_neon = _load("Neon").trim_axes(method="x-axis", boundaries=(100, 3500)) + papp = _build_measurement_papp( + [(spe_neon.x, spe_neon.y, "RAW_DATA", "reference_neon", "cm-1")], + meta=None, instrument=("TestCo", "ModelX"), wavelength=LASER_WL) + assert papp is not None + + nx_root = nx.NXroot() + papp.to_nexus(nx_root) + entry = next(iter(nx_root.entries.values())) + assert "instrument" in entry + path = str(tmp_path / "papp_only.nxs") + nx_root.save(path, mode="w") + + import h5py + with h5py.File(path, "r") as f: + groups = [] + f.visititems(lambda name, obj: groups.append(name)) + assert any("reference_neon" in g for g in groups), ( + f"reference_neon effect not found anywhere in the written tree: {groups}") + + +def _entry_path(h5file): + """The one real NXentry in the file. pyambit's to_nexus writes it at a leading-slash + path derived from (provider, papp.uuid), not a fixed "/entry" -- discover it instead + of assuming the name (see test_build_measurement_papp_entry_is_discoverable).""" + names = [k for k, v in h5file.items() + if v.attrs.get("NX_class") in (b"NXentry", "NXentry")] + assert len(names) == 1, f"expected exactly one NXentry, found {names}" + return names[0] + + +def test_nexusformat_can_open(calmodel, tmp_path): + """Structural sanity: the file is a valid NeXus tree, not just valid HDF5. + + nexusformat consumes the ``NX_class`` HDF5 attribute into its own ``.nxclass`` + property (it is not left visible via ``.attrs``), so a passing ``.nxclass`` check here + is what proves nexusformat parsed our groups as real NeXus classes rather than opaque + HDF5 groups. + """ + nexusformat = pytest.importorskip("nexusformat.nexus") + path = _write(calmodel, tmp_path) + root = nexusformat.nxload(path) + entry = next(e for e in root.entries.values() if e.nxclass == "NXentry") + assert entry.instrument.nxclass == "NXinstrument" + assert entry.instrument.calibration_x_0.nxclass == "NXcalibration" + assert entry.calibration_curve.nxclass == "NXdata" + + +def test_calibration_object_json_roundtrips(calmodel, spe_pst, tmp_path): + """The full model, reconstructed from the file alone, reproduces the calibrated axis.""" + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + doc = json.loads(f[f"{entry}/calibration_model/data"][()]) + reconstructed = CalibrationModel.from_dict(doc["model"]) + + orig = calmodel.apply_calibration_x(spe_pst, spe_units="cm-1").x + rt = reconstructed.apply_calibration_x(spe_pst, spe_units="cm-1").x + np.testing.assert_allclose(rt, orig, rtol=0, atol=1e-9) + + +def test_calibration_parameters_written(calmodel, tmp_path): + """The reconstructable fit parameters are real numeric datasets, not just JSON text.""" + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + group = f[f"{entry}/instrument/calibration_x_0/calibration_parameters"] + kind = calmodel.components[0].to_dict()["model"]["type"] + if kind == "CustomPolyInterpolator": + assert "coef" in group and group["coef"].shape[0] >= 1 + assert int(group["degree"][()]) >= 1 + assert group["x_min"][()] < group["x_max"][()] + elif kind == "CustomPChipInterpolator": + assert group["knots_x"].shape == group["knots_y"].shape + assert group["knots_x"].shape[0] > 1 + else: + pytest.fail(f"unexpected interpolator kind {kind!r}") + + si_group = f[f"{entry}/instrument/calibration_x_1/calibration_parameters"] + assert si_group["si_peak_nm"][()] > 0 + + +def test_original_and_calibrated_axis_same_length(calmodel, tmp_path): + """The sampled curve datasets satisfy the NXDL `ncal` symbol constraint (equal length).""" + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + original = f[f"{entry}/instrument/calibration_x_0/original_axis"] + calibrated = f[f"{entry}/instrument/calibration_x_0/calibrated_axis"] + assert original.shape == calibrated.shape + assert original.shape[0] == 101 # npoints passed to _write + + +def test_anchors_not_in_calibrated_axis(calmodel, tmp_path): + """Matched-peak anchors (a different length than the sampled curve) live in their own + group, not mixed into original_axis/calibrated_axis.""" + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + cal = f[f"{entry}/instrument/calibration_x_0"] + assert "anchors" in cal + anchors = cal["anchors"] + assert set(anchors.keys()) >= {"measured", "reference", "inlier"} + assert anchors["measured"].shape == anchors["reference"].shape + + +def test_reference_spectra_written(calmodel, spe_neon, spe_sil, tmp_path): + """Regression guard: the calibrant x AND y arrays are actually persisted (the failure + mode this guards against is ramanchada2.io.HSDS.write_nexus silently dropping the axis + via require_group(..., data=x) — a group has no data= kwarg). Written through pyambit's + spe2ambit/process_pa, so the effect lands under a RAW_DATA group, named by nx_name.""" + path = _write(calmodel, tmp_path, spe_neon=spe_neon, spe_neon_units="cm-1", + spe_silicon=spe_sil, spe_silicon_units="cm-1") + with h5py.File(path, "r") as f: + entry = _entry_path(f) + neon = f[f"{entry}/RAW_DATA/reference_neon_1"] + np.testing.assert_allclose(neon["raman_shift"][()], np.asarray(spe_neon.x)) + np.testing.assert_allclose(neon["y"][()], np.asarray(spe_neon.y)) + + sil = f[f"{entry}/RAW_DATA/reference_silicon_2"] + np.testing.assert_allclose(sil["raman_shift"][()], np.asarray(spe_sil.x)) + np.testing.assert_allclose(sil["y"][()], np.asarray(spe_sil.y)) + + +def test_reference_spectra_omitted_when_not_given(calmodel, tmp_path): + """No invented data: no RAW_DATA reference-spectrum effect exists if the caller didn't + supply calibrant spectra (the x-calibration curve itself still gets its own + non-RAW_DATA group, so RAW_DATA is absent entirely rather than merely smaller).""" + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + assert "RAW_DATA" not in f[entry] + + +def test_ycal_component_written_when_given(calmodel, ycal, spe_pst, tmp_path): + path = _write(calmodel, tmp_path, ycal_component=ycal) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + ycal_group = f[f"{entry}/instrument/calibration_y"] + assert ycal_group.attrs["physical_quantity"] == "relative intensity" + doc = json.loads(ycal_group["calibration_object/data"][()]) + reconstructed = YCalibrationComponent.from_dict(doc) + spe = spe_pst.trim_axes(method="x-axis", boundaries=ycal.ref.raman_shift) + np.testing.assert_allclose( + reconstructed.process(spe).y, ycal.process(spe).y, rtol=0, atol=1e-9) + + +def test_ycal_plottable_curves_present(calmodel, ycal, tmp_path): + """Regression: the y-calibration must be visually discoverable, not just present as + fit parameters/JSON. The x side already had a plottable calibration_curve NXdata; + the y side previously had none, which is why a user opening the file couldn't find + the intensity calibration at all despite the model being fully written.""" + path = _write(calmodel, tmp_path, ycal_component=ycal) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + assert "calibration_curve_y" in f[entry], ( + "no plottable NXdata for the y (relative-intensity) calibration curve") + y_curve = f[f"{entry}/calibration_curve_y"] + assert y_curve.attrs["NX_class"] == "NXdata" + assert y_curve.attrs["signal"] == "intensity_factor" + assert y_curve["intensity_factor"].shape[0] > 1 + + # the measured SRM response and the certificate's analytic response are both + # written as their own RAW_DATA-adjacent effects (via spe2ambit), not just + # embedded numbers in calibration_object's JSON + raw_names = list(f[f"{entry}/Y_CALIBRATION_MEASURED"].keys()) \ + if "Y_CALIBRATION_MEASURED" in f[entry] else [] + cert_names = list(f[f"{entry}/Y_CALIBRATION_CERTIFICATE"].keys()) \ + if "Y_CALIBRATION_CERTIFICATE" in f[entry] else [] + assert raw_names, "measured SRM response effect not written" + assert cert_names, "certificate response effect not written" + + +def test_ycal_component_omitted_when_not_given(calmodel, tmp_path): + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + assert "calibration_y" not in f[f"{entry}/instrument"] + + +def test_default_attr_points_at_calibration_curve(calmodel, tmp_path): + path = _write(calmodel, tmp_path) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + default_name = f[entry].attrs["default"] + assert default_name in f[entry] + assert f[entry][default_name].attrs["NX_class"] == "NXdata" + + +def test_instrument_metadata_written_and_nan_dropped(calmodel, tmp_path): + """instrument_make/instrument_model become pyambit's typed (vendor, model) device + identity; other instrument keys (grating) route through configure_papp's backward- + compat table into typed fields or the generic parameters bucket. NaN values must not + reach either -- they're dropped before crossing into pyambit's meta dict.""" + path = _write(calmodel, tmp_path, instrument={ + "instrument_make": "TestCo", "instrument_model": "ModelX", + "grating": float("nan"), "laser_wl": 785, + }) + with h5py.File(path, "r") as f: + entry = _entry_path(f) + device = f[f"{entry}/instrument/device_information"] + assert device["vendor"][()].decode() == "TestCo" + assert device["model"][()].decode() == "ModelX" + # grating was NaN -> dropped; must not appear anywhere as a written value + params = f[f"{entry}/parameters"] if "parameters" in f[entry] else {} + assert "grating" not in params From 4aa3f4c68861ba21be7e38f4f257bf0d382990ef Mon Sep 17 00:00:00 2001 From: vedina Date: Mon, 5 Oct 2026 22:17:03 +0300 Subject: [PATCH 2/2] poetry lock updated --- poetry.lock | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/poetry.lock b/poetry.lock index eb1925a..8f0d444 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1,4 +1,4 @@ -# This file is automatically @generated by Poetry 2.4.1 and should not be changed by hand. +# This file is automatically @generated by Poetry 2.2.1 and should not be changed by hand. [[package]] name = "annotated-types" @@ -3340,4 +3340,4 @@ type = ["pytest-mypy (>=1.0.1) ; platform_python_implementation != \"PyPy\""] [metadata] lock-version = "2.1" python-versions = ">=3.10,<3.14" -content-hash = "cfce4d7c2f404ae0f70a59b8e42989beffcbe4be0cd2296e27da42d551cf4f68" +content-hash = "b9f0c31af3385c17d4ed47c437567cb07d25509dc19a4b553e694264be26a195"