diff --git a/.github/workflows/build-and-publish-docs-on-dispatch.yml b/.github/workflows/build-and-publish-docs-on-dispatch.yml index a60668b..1be1b4f 100644 --- a/.github/workflows/build-and-publish-docs-on-dispatch.yml +++ b/.github/workflows/build-and-publish-docs-on-dispatch.yml @@ -13,6 +13,7 @@ jobs: uses: scikit-package/release-scripts/.github/workflows/_release-docs.yml@v0 with: project: pyobjcryst + submodules: "true" c_extension: true headless: false python_version: ${{ fromJSON(needs.get-python-version.outputs.latest_python_version) }} diff --git a/.github/workflows/build-wheel-release-upload.yml b/.github/workflows/build-wheel-release-upload.yml index 396b9ec..42b91f7 100644 --- a/.github/workflows/build-wheel-release-upload.yml +++ b/.github/workflows/build-wheel-release-upload.yml @@ -27,6 +27,7 @@ jobs: uses: scikit-package/release-scripts/.github/workflows/_build-wheel-release-upload.yml@v0 with: project: pyobjcryst + submodules: "true" c_extension: true maintainer_github_username: sbillinge secrets: diff --git a/.github/workflows/matrix-and-codecov.yml b/.github/workflows/matrix-and-codecov.yml index 9a3b8e2..f02cc91 100644 --- a/.github/workflows/matrix-and-codecov.yml +++ b/.github/workflows/matrix-and-codecov.yml @@ -15,6 +15,7 @@ jobs: uses: scikit-package/release-scripts/.github/workflows/_matrix-and-codecov-on-merge-to-main.yml@v0 with: project: pyobjcryst + submodules: "true" c_extension: true headless: false secrets: diff --git a/.github/workflows/tests-on-pr.yml b/.github/workflows/tests-on-pr.yml index 6293743..2480080 100644 --- a/.github/workflows/tests-on-pr.yml +++ b/.github/workflows/tests-on-pr.yml @@ -9,6 +9,7 @@ jobs: uses: scikit-package/release-scripts/.github/workflows/_tests-on-pr.yml@v0 with: project: pyobjcryst + submodules: "true" c_extension: true headless: false secrets: diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b60917a..beb8815 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,3 +1,5 @@ +# Vendored third-party sources (objcryst submodule, cctbx, newmat) +exclude: "^extern/" default_language_version: python: python3 ci: diff --git a/nanobind_migration_notes.md b/nanobind_migration_notes.md index ffc2e35..57eb662 100644 --- a/nanobind_migration_notes.md +++ b/nanobind_migration_notes.md @@ -6,6 +6,7 @@ gaps listed below — it is not "fully done" in the sense of zero open items.** What's done and verified: + - The C++ binding layer is fully ported from Boost.Python + SCons to nanobind + CMake + scikit-build-core. `src/pyobjcryst/*.py` (the pure-Python layer) and `tests/*.py` are untouched — nothing there is Boost.Python-specific. @@ -21,6 +22,7 @@ What's done and verified: - **161/161 tests pass** (`pytest tests/ -q`). Known, accepted gaps — real, not oversights, but still open: + - **Python-side method override does not work across the `RefinableObj` virtual-inheritance seam.** A Python subclass of `Atom`, `Molecule`, `ScatteringPowerAtom`, `DiffractionDataSingleCrystal`, `PowderPatternBackground`, etc. that overrides an inherited @@ -28,7 +30,7 @@ Known, accepted gaps — real, not oversights, but still open: C++-internal virtual dispatch — only plain forwarding calls were restored, not overridability. This was an explicit scope decision (general override support was confirmed not needed), not an oversight, but it is a real regression vs. Boost.Python's `bases`, which handled - this transparently. Override support that *does* still work: `RefinableObj`'s own ~15 virtuals + this transparently. Override support that _does_ still work: `RefinableObj`'s own ~15 virtuals when subclassing `RefinableObj` directly; `Scatterer`/`ScatteringPower`/`Restraint`'s own (non-inherited) virtuals; and `UpdateDisplay` specifically on `Crystal`, `PowderPattern`, and `MonteCarloObj` (all three reach `RefinableObj` through a non-virtual chain, so no seam exists @@ -153,18 +155,19 @@ after every test run — pre-existing, tracked in [Status](#status-2026-08-12) a ### Build troubleshooting -| Symptom | Cause | Fix | -|---|---|---| -| `nanobind not found` | not installed | `conda install nanobind` | -| `ObjCryst submodule empty` | not initialized | `git submodule update --init --recursive` | -| `version.h: No such file` | CMake config not run | `pip install -e . --no-build-isolation` | -| Segfault / bus error in a test | pointer-offset bug (virtual inheritance) | check whether the crashing class's C++ header declares `virtual` inheritance anywhere in its ancestry; see [Virtual Inheritance](#virtual-inheritance) | -| `RuntimeError: std::bad_cast`, no useful location | usually `nb::cast` on a `bytes` object | check every `.attr("read")()` call site — see [Python file-like input](#python-file-like-input-str-vs-bytes) | -| Binary needs a newer glibc than the host ships (`GLIBC_ABI_GNU2_TLS`) | conda-forge GCC's `gnu2` TLS dialect | build with the system compiler: `CC=/usr/bin/gcc CXX=/usr/bin/g++ pip install -e . --no-build-isolation` | -| gdb: `(no debugging symbols found)` | scikit-build-core strips symbols on install | `pip install -e . --no-build-isolation -Ccmake.define.CMAKE_BUILD_TYPE=RelWithDebInfo`, then copy the unstripped `.so` over the installed one | -| Values off in the last few digits | `REAL` (`ObjCryst/ObjCryst/General.h`) defaults to `float`; must be forced to `double` | confirm `CMakeLists.txt` sets `target_compile_definitions(ObjCryst PUBLIC REAL=double ...)`, and that any helper uses the `CrystVector_REAL`/`CrystMatrix_REAL` macros rather than a hardcoded ``/`` | +| Symptom | Cause | Fix | +| --------------------------------------------------------------------- | -------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `nanobind not found` | not installed | `conda install nanobind` | +| `ObjCryst submodule empty` | not initialized | `git submodule update --init --recursive` | +| `version.h: No such file` | CMake config not run | `pip install -e . --no-build-isolation` | +| Segfault / bus error in a test | pointer-offset bug (virtual inheritance) | check whether the crashing class's C++ header declares `virtual` inheritance anywhere in its ancestry; see [Virtual Inheritance](#virtual-inheritance) | +| `RuntimeError: std::bad_cast`, no useful location | usually `nb::cast` on a `bytes` object | check every `.attr("read")()` call site — see [Python file-like input](#python-file-like-input-str-vs-bytes) | +| Binary needs a newer glibc than the host ships (`GLIBC_ABI_GNU2_TLS`) | conda-forge GCC's `gnu2` TLS dialect | build with the system compiler: `CC=/usr/bin/gcc CXX=/usr/bin/g++ pip install -e . --no-build-isolation` | +| gdb: `(no debugging symbols found)` | scikit-build-core strips symbols on install | `pip install -e . --no-build-isolation -Ccmake.define.CMAKE_BUILD_TYPE=RelWithDebInfo`, then copy the unstripped `.so` over the installed one | +| Values off in the last few digits | `REAL` (`ObjCryst/ObjCryst/General.h`) defaults to `float`; must be forced to `double` | confirm `CMakeLists.txt` sets `target_compile_definitions(ObjCryst PUBLIC REAL=double ...)`, and that any helper uses the `CrystVector_REAL`/`CrystMatrix_REAL` macros rather than a hardcoded ``/`` | Useful gdb recipe for a segfault/bad_cast inside the extension: + ```bash gdb -q --args python -m pytest tests/test_whatever.py -k test_name (gdb) catch throw # for std::bad_cast / C++ exceptions @@ -182,7 +185,7 @@ or extending any class near `RefinableObj`, `Scatterer`, `ScatteringPower`, `Sca ### The core problem -ObjCryst++ uses *virtual* inheritance in several places: `Scatterer : virtual public +ObjCryst++ uses _virtual_ inheritance in several places: `Scatterer : virtual public RefinableObj`, `ScatteringPower : virtual public RefinableObj`, `ScatteringData : virtual public RefinableObj`, `PowderPatternComponent : virtual public RefinableObj`, and a second hop for `ScatteringPowerAtom : virtual public ScatteringPower` and `GlobalScatteringPower : virtual public @@ -190,8 +193,8 @@ ScatteringPower`. nanobind's automatic base↔derived pointer adjustment (used by `nb::class_(...)`) assumes a **fixed, statically computable** offset between `Derived*` and `Base*`. That's only -valid for *non-virtual* inheritance. Declaring the nanobind base across a virtual link does **not** -raise an error — it silently computes a *wrong* offset, which can look fine for fields at offset 0 +valid for _non-virtual_ inheritance. Declaring the nanobind base across a virtual link does **not** +raise an error — it silently computes a _wrong_ offset, which can look fine for fields at offset 0 and then segfault or return garbage the moment a method touching state further into the object is called. @@ -203,6 +206,12 @@ is fine and should still use nanobind's automatic base declaration (e.g. `Scatte public ScatteringPower`, no `virtual`, is safely `nb::class_`). +The same applies to _non-virtual multiple inheritance_ where the base is not the first one, e.g. +`TextureMarchDollase : public ScatteringCorr, public RefinableObj`: `RefinableObj` sits at a +non-zero offset, and `nb::class_` silently returns garbage +(`GetNbPar()` was 0). Register it standalone, call `bind_refinableobj_forwarding()`, and add it to +`extractRefinableObjArg`. + ### Four ways this bites, and the fix for each **1. Return direction** — a function returns a base-typed reference/pointer @@ -218,14 +227,15 @@ nb::object wrapScatteringPowerReturn(const ScatteringPower* sp, nb::handle paren return nb::cast(const_cast(sp), nb::rv_policy::reference_internal, parent); } ``` -Applies even when the returned type is an *intermediate* base, not the ultimate leaf (e.g. + +Applies even when the returned type is an _intermediate_ base, not the ultimate leaf (e.g. `PowderPattern::GetPowderPatternComponent(i)` returns `PowderPatternComponent&`, whose dynamic type is `PowderPatternDiffraction`) — every function returning a shared ancestor type anywhere in a virtual hierarchy needs this, not just the one function returning the immediate leaf type. -**2. Argument direction** — a virtually-derived leaf is registered *without* a declared nanobind +**2. Argument direction** — a virtually-derived leaf is registered _without_ a declared nanobind base (per the rule above), so nanobind refuses it wherever a C++ function's parameter type is the -*base* (`void AddAtom(..., const ScatteringPower* pow, ...)` rejects a Python `ScatteringPowerAtom` +_base_ (`void AddAtom(..., const ScatteringPower* pow, ...)` rejects a Python `ScatteringPowerAtom` with `TypeError: incompatible function arguments`). Fix: accept `nb::object` and resolve manually, trying the concrete leaf type first: @@ -238,13 +248,13 @@ const ScatteringPower* extractScatteringPowerArg(nb::object obj) { ``` **3. Generic "any `RefinableObj`" argument functions** — the same argument-direction problem, but -for functions generic over *any* `RefinableObj` (`MonteCarloObj::AddRefinableObj`, +for functions generic over _any_ `RefinableObj` (`MonteCarloObj::AddRefinableObj`, `LSQNumObj::SetRefinedObj`, `RefinableObj::AddPar`): nanobind has no conversion path from a Python `DiffractionDataSingleCrystal`/`Atom`/`Molecule`/etc. to `RefinableObj&`, since none of those concrete classes declare that base. Fix: a shared helper (`extractRefinableObjArg`, in `helpers_nb.hpp`/`nb_refinableobj.cpp`) trying the safe cast first, then falling back through every known virtually-linked type with `static_cast` (always well-defined at compile time — the problem -is only ever nanobind's own *runtime* offset table, never a real `static_cast`): +is only ever nanobind's own _runtime_ offset table, never a real `static_cast`): ```cpp RefinableObj& extractRefinableObjArg(nb::object obj) { @@ -258,11 +268,12 @@ RefinableObj& extractRefinableObjArg(nb::object obj) { } ``` -**4. Under-declaring a *safe* base** — the mirror-image mistake: forgetting to declare a nanobind +**4. Under-declaring a _safe_ base** — the mirror-image mistake: forgetting to declare a nanobind base across a link that's perfectly safe (plain, non-virtual). Symptom looks identical -(`TypeError: incompatible function arguments`) but the cause is a *missing* declaration, not a +(`TypeError: incompatible function arguments`) but the cause is a _missing_ declaration, not a dangerous one. Fix: register the base type too (it needs no bindings of its own, just an internal name), then declare the relationship normally: + ```cpp nb::class_(m, "_PeakListBase"); // internal-only nb::class_(m, "PeakList"); // now declared @@ -281,28 +292,29 @@ Fix: two reusable template helpers in `helpers_nb.hpp`, `bind_refinableobj_forwa `bind_scatteringpower_forwarding()`, each binding one method-set's full surface via plain forwarding `.def()` calls (`&T::Method`, or a lambda for anything needing a helper/numpy conversion). This is safe despite the virtual base because `T` is fully known at compile time at -the call site, so the *compiler* resolves the pointer adjustment — exactly the mechanism +the call site, so the _compiler_ resolves the pointer adjustment — exactly the mechanism Boost.Python relied on, nothing to do with nanobind's own (virtual-base-unsafe) casting. Called once each in `nb_scatterer.cpp`, `nb_scatteringpower.cpp`, `nb_scatteringdata.cpp`, and `nb_powderpatterncomponent.cpp` (cascades automatically to every non-virtually-inherited child — `Atom`, `Molecule`, `ZScatterer`, `PowderPatternBackground`, `PowderPatternDiffraction`, etc. — via ordinary Python MRO); `ScatteringPowerAtom` and `GlobalScatteringPower` need both helpers called -directly on them, since they have a *second* virtual-inheritance hop of their own. +directly on them, since they have a _second_ virtual-inheritance hop of their own. **This restores callability, not overridability** — see [Status](#status-2026-08-12) for why that distinction matters and is a deliberate, accepted scope limit (nanobind's `NB_TRAMPOLINE` macro -requires its argument to be a *direct* base, a hard C++ rule that blocks chained/templated +requires its argument to be a _direct_ base, a hard C++ rule that blocks chained/templated trampolines the way pybind11 supports them — confirmed by compiler error in a standalone repro, not a nanobind design choice that could be worked around). Two real bugs turned up while adding test coverage for this fix, both worth remembering as general lessons: + - A binding that reaches a protected member via a trampoline/alias helper class (here, `PyScatterer::GetClockScattCompList()`, added to reach `Scatterer`'s protected `mClockScattCompList`) must take the **real base type** (`Scatterer&`) as its bound parameter, not the trampoline/alias type (`PyScatterer&`) — a plain `Atom`/`Molecule` instance is never actually constructed as the trampoline, so binding against the alias type raises `TypeError` for - every ordinary instance. Cast through the alias *inside* the lambda instead: + every ordinary instance. Cast through the alias _inside_ the lambda instead: ```cpp .def("GetClockScattCompList", [](Scatterer& s) -> const RefinableObjClock& { @@ -310,7 +322,7 @@ lessons: }, nb::rv_policy::reference_internal) ``` - `bind_scatteringpower_forwarding()` forwarded `GetBij`/`SetBij` but missed the `B11`/`B12`/ - `B13`/`B22`/`B23`/`B33`/`Biso` *properties* that `nb_scatteringpower.cpp` defines separately via + `B13`/`B22`/`B23`/`B33`/`Biso` _properties_ that `nb_scatteringpower.cpp` defines separately via file-local templates — a reminder that a forwarding helper needs to be checked against the full `dir()` of the canonical class, not just its method list, when a class also exposes properties. @@ -321,9 +333,11 @@ lessons: ### `None` argument handling nanobind rejects `None` by default for any typed `nb::arg`. Allow it explicitly: + ```cpp .def("AddScatterer", &Crystal::AddScatterer, nb::arg("pScatt").none(), nb::arg("ownsScatt") = true) ``` + `.none()` is required even for generic `nb::object` parameters — it isn't automatic just because the C++ parameter type is generic. @@ -332,6 +346,7 @@ the C++ parameter type is generic. `REAL* derivpar` in `RecUnitCell::hkl2d`) — nanobind's arithmetic `type_caster` rejects `None` outright before `.none()` on the outer `nb::arg` ever gets a chance to matter. Fix: wrap in a helper taking `nb::object` and bridge manually: + ```cpp float _hkl2d(const RecUnitCell& r, float h, float k, float l, nb::object derivpar, unsigned int derivhkl) { if (derivpar.is_none()) return r.hkl2d(h, k, l, nullptr, derivhkl); @@ -351,8 +366,8 @@ produce wrong-precision code. ### Lambda return type silently deduced by value A `.def(...)` that wraps a reference-returning helper in a lambda **without an explicit trailing -return type** deduces `auto`, which strips the reference and returns *by value* — silently -*copying* the object, even with `nb::rv_policy::reference_internal` set (that policy only controls +return type** deduces `auto`, which strips the reference and returns _by value_ — silently +_copying_ the object, even with `nb::rv_policy::reference_internal` set (that policy only controls how nanobind wraps whatever the lambda actually returns; it can't recover a reference from an already-copied value). This caused a real bug: `Molecule.AddRigidGroup`'s lambda copied the `RigidGroup` it had just added, so `Molecule.RemoveRigidGroup` (which finds the group to remove by @@ -371,18 +386,18 @@ confirmed by a plain, non-pytest script that creates a `Molecule`, calls `AddBon all references, and calls `gc.collect()` twice: still leaks 6 instances at exit. Root cause: nanobind types don't participate in Python's cyclic GC by default (no `tp_traverse`/`tp_clear`), so a bidirectional `keep_alive` between two live objects — parent keeps child alive (a `.def()` -declares `keep_alive<1,N>()`) *and* child keeps parent alive (a `.def()` returns via +declares `keep_alive<1,N>()`) _and_ child keeps parent alive (a `.def()` returns via `rv_policy::reference_internal`, nanobind's shorthand for `keep_alive<0,1>()`) — forms a cycle nothing can ever collect. Fixed for `Molecule`↔`MolAtom`: removed the redundant Molecule-keeps-atom direction from `AddBond`/`AddBondAngle`/`AddDihedralAngle` (safe because `MolAtom` has no standalone Python-owned existence — its only constructor is a copy constructor, so it's always reached via -`GetAtom()`/`AddAtom()`, which already protects it via the *other* direction). +`GetAtom()`/`AddAtom()`, which already protects it via the _other_ direction). **Not safe to fix the same way**, and still open: `Crystal`↔`Scatterer`/`Molecule`/ `ScatteringPower` (and possibly `RefinableObj`↔`RefinablePar`). Unlike `MolAtom`, these classes -*do* have real standalone Python-owned constructors (`sp = ScatteringPowerAtom("C", "C")` with no +_do_ have real standalone Python-owned constructors (`sp = ScatteringPowerAtom("C", "C")` with no parent) — the forward direction is protecting a real gap (between construction and being added to a `Crystal`, nothing else keeps the child's wrapper alive), so removing it would trade a leak for a use-after-free. The real fix is nanobind's `tp_traverse`/`tp_clear` mechanism, which likely requires @@ -396,6 +411,7 @@ the parent — that combination is the signature to look for. Any binding calling `.attr("read")()` on a caller-supplied Python object and casting straight to `std::string` breaks the moment the file was opened in binary mode (`.read()` returns `bytes`, `nb::cast` throws `std::bad_cast`). Use a shared helper: + ```cpp inline std::string read_pyfile_to_string(nb::object input) { nb::object data = input.attr("read")(); @@ -437,6 +453,7 @@ this fix unless it's re-applied or merged upstream. ## Binding Patterns Quick Reference **Simple class, no virtuals:** + ```cpp nb::class_(m, "MyClass") .def(nb::init<>()) @@ -445,8 +462,9 @@ nb::class_(m, "MyClass") ; ``` -**Class with overridable virtuals (trampoline) — only where `Base` is a *direct*, non-virtual +**Class with overridable virtuals (trampoline) — only where `Base` is a _direct_, non-virtual base of the trampoline, per [Virtual Inheritance](#virtual-inheritance):** + ```cpp struct PyMyClass : MyClass { NB_TRAMPOLINE(MyClass, 2); // 2 virtual methods @@ -457,6 +475,7 @@ nb::class_(m, "MyClass").def(nb::init<>()); ``` **String representation** (many ObjCryst objects have a useful `operator<<`): + ```cpp .def("__str__", [](const T& obj){ return obj_str(obj); }) // obj_str() in helpers_nb.hpp, mutes debug spam .def("__repr__", [](const T& obj){ return obj_str(obj); }) @@ -464,6 +483,7 @@ nb::class_(m, "MyClass").def(nb::init<>()); **Method overload resolution** — nanobind can't distinguish C++ overloads at binding time; use a lambda or a pointer-to-member-function cast: + ```cpp .def("SetX", nb::overload_cast(&T::SetX)) ``` @@ -476,13 +496,13 @@ handling, container conversion, lambda return types, `keep_alive`, and file-like ## Architecture Decisions -| Decision | Rationale | -|---|---| -| Single `_pyobjcryst` extension | simpler than per-class extensions; nanobind modules can be large | -| Static linking (newmat, cctbx, ObjCryst) | avoids runtime library conflicts on Windows/macOS wheels | -| scikit-build-core + CMake | modern standard for C++ Python extensions, replaces SCons | -| ObjCryst++ as a submodule | tracks upstream, simpler than a fork | -| Vendored newmat/cctbx (not conda-forge `cctbx-base`) | see below | +| Decision | Rationale | +| ---------------------------------------------------- | ---------------------------------------------------------------- | +| Single `_pyobjcryst` extension | simpler than per-class extensions; nanobind modules can be large | +| Static linking (newmat, cctbx, ObjCryst) | avoids runtime library conflicts on Windows/macOS wheels | +| scikit-build-core + CMake | modern standard for C++ Python extensions, replaces SCons | +| ObjCryst++ as a submodule | tracks upstream, simpler than a fork | +| Vendored newmat/cctbx (not conda-forge `cctbx-base`) | see below | **Vendored `cctbx` vs. conda-forge's `cctbx-base`** (investigated, not adopted): `cctbx-base` does contain the right modules, compiled and working, but isn't structured for third-party C++ @@ -499,17 +519,19 @@ depending on `cctbx-base` as it currently ships. ## Pitfalls Checklist **Before writing/extending a binding for class `X`:** -- [ ] Grep `X`'s header and its whole ancestry (not just the direct parent) for `: *virtual - public` — if found, do **not** declare that base in `nb::class_`; use the manual - `dynamic_cast`/`nb::cast` patterns instead + +- [ ] Grep `X`'s header and its whole ancestry (not just the direct parent) for + `: *virtual public` — if found, do **not** declare that base in + `nb::class_`; use the manual `dynamic_cast`/`nb::cast` patterns instead **While writing:** + - [ ] Every optional pointer/object argument needs `.none()` on its `nb::arg(...)`, including plain `nb::object` params — except pointers to a plain arithmetic type, where `.none()` has no effect; wrap those in an `nb::object`-taking helper instead - [ ] Every function returning a base-typed reference/pointer whose dynamic type is commonly a virtually-derived leaf needs the `dynamic_cast`-and-recast treatment -- [ ] Every function generically *accepting* a base type as an argument needs its argument routed +- [ ] Every function generically _accepting_ a base type as an argument needs its argument routed through a helper like `extractRefinableObjArg()` - [ ] Any `.def(...)` combining a lambda with `nb::rv_policy::reference`/`reference_internal` needs an explicit `-> T&`/`-> T*` trailing return type on the lambda @@ -518,10 +540,11 @@ depending on `cctbx-base` as it currently ships. - [ ] Any container type (`CrystVector`, `CrystMatrix`) crossing the boundary should use the `CrystVector_REAL`/`CrystMatrix_REAL` macros, never a hardcoded ``/`` - [ ] A binding reaching a protected/private member via a trampoline/alias helper class must take - the *real* base type as its bound parameter, not the trampoline/alias type — otherwise it + the _real_ base type as its bound parameter, not the trampoline/alias type — otherwise it raises `TypeError` on every instance not literally constructed as that alias **After writing, before declaring done:** + - [ ] Rebuilt and re-ran the specific test file, then the full suite, and confirmed no new failures - [ ] If the class exposes properties as well as methods, checked the forwarding/coverage against the full `dir()` of the canonical class, not just its method list diff --git a/news/march-dollase-texture.rst b/news/march-dollase-texture.rst new file mode 100644 index 0000000..3a891a0 --- /dev/null +++ b/news/march-dollase-texture.rst @@ -0,0 +1,23 @@ +**Added:** + +* Expose the March-Dollase preferred-orientation correction TextureMarchDollase (AddPhase, SetPhasePar, DeletePhase and getters) and PowderPatternDiffraction.GetTextureMarchDollase. SetPhasePar and DeletePhase rely on the upstream ObjCryst++ fix (vincefn/objcryst#105) which recalculates the powder pattern when the texture changes. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/news/tof-profile.rst b/news/tof-profile.rst new file mode 100644 index 0000000..ad9d2c6 --- /dev/null +++ b/news/tof-profile.rst @@ -0,0 +1,23 @@ +**Added:** + +* Expose the time-of-flight reflection profile ReflectionProfileDoubleExponentialPseudoVoigt. + +**Changed:** + +* PowderPatternDiffraction.SetProfile gives a ReflectionProfileDoubleExponentialPseudoVoigt the unit cell of the component's crystal. + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/extensions/nb_indexing.cpp b/src/extensions/nb_indexing.cpp index e06ace5..44bf2a3 100644 --- a/src/extensions/nb_indexing.cpp +++ b/src/extensions/nb_indexing.cpp @@ -30,7 +30,7 @@ namespace { // with "incompatible function arguments", even though None is exactly what a // Python caller must pass (a Python float is immutable, so there is no way for // the derivative to be written back to it regardless of binding technology -- -// this matches the pre-existing Boost.Python behavor, where derivpar was only +// this matches the pre-existing Boost.Python behavior, where derivpar was only // ever meaningfully used from within C++ itself, e.g. LSQ fitting internals). // Fixed by taking derivpar as nb::object and manually bridging None -> nullptr. float _hkl2d(const RecUnitCell& r, const float h, const float k, const float l, diff --git a/src/extensions/nb_powderpatterndiffraction.cpp b/src/extensions/nb_powderpatterndiffraction.cpp index 92df37c..df8627f 100644 --- a/src/extensions/nb_powderpatterndiffraction.cpp +++ b/src/extensions/nb_powderpatterndiffraction.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include "helpers_nb.hpp" @@ -18,6 +19,38 @@ namespace nb = nanobind; using namespace ObjCryst; +namespace { + +void SetProfileCopy(PowderPatternDiffraction& pdiff, + const ReflectionProfile& profile) +{ + ReflectionProfile* copy = profile.CreateCopy(); + // The TOF profile needs a unit cell to compute d_hkl, as in + // PowderPatternDiffraction::SetCrystal and XMLInput. + auto* tofprofile = + dynamic_cast(copy); + if (tofprofile != nullptr && pdiff.HasCrystal()) + tofprofile->SetUnitCell(pdiff.GetCrystal()); + pdiff.SetProfile(copy); +} + +// PowderPatternDiffraction has no public accessor for its texture +// correction, but registers it as a sub-object. +TextureMarchDollase& GetTextureMarchDollase(PowderPatternDiffraction& pdiff) +{ + ObjRegistry& subobjs = pdiff.GetSubObjRegistry(); + for (long i = 0; i < subobjs.GetNb(); ++i) + { + auto* texture = dynamic_cast(&subobjs.GetObj(i)); + if (texture != nullptr) return *texture; + } + throw ObjCrystException("PowderPatternDiffraction has no " + "TextureMarchDollase sub-object. Please report " + "this as a pyobjcryst bug."); +} + +} // namespace + void wrap_powderpatterndiffraction(nb::module_& m) { nb::enum_(m, "ReflectionProfileType") @@ -72,9 +105,7 @@ void wrap_powderpatterndiffraction(nb::module_& m) nb::overload_cast<>(&PowderPatternDiffraction::GetProfile), nb::rv_policy::reference_internal) .def("SetProfile", - [](PowderPatternDiffraction& d, ReflectionProfile& p) { - d.SetProfile(p.CreateCopy()); - }, + &SetProfileCopy, nb::arg("profile"), "Install an independent copy of the given profile.") .def("SetExtractionMode", @@ -94,5 +125,31 @@ void wrap_powderpatterndiffraction(nb::module_& m) [](PowderPatternDiffraction& d){ return crystvec_to_array(d.GetFhklObsSq()); }) .def("X2XCorrPhase", &PowderPatternDiffraction::X2XCorrPhase, "Apply the flat-detector displacement correction for this phase's offset.") + .def("GetTextureMarchDollase", &GetTextureMarchDollase, + nb::rv_policy::reference_internal, + "Return the March-Dollase preferred-orientation correction.") + ; + + // TextureMarchDollase : public ScatteringCorr, public RefinableObj -- + // RefinableObj is a second base at a non-zero offset, which nanobind's + // base casting cannot handle (same failure mode as the virtual bases in + // nanobind_migration_notes.md). Register it standalone and forward the + // RefinableObj surface through T instead. + auto texture = nb::class_(m, "TextureMarchDollase"); + bind_refinableobj_forwarding(texture); + texture + .def("AddPhase", &TextureMarchDollase::AddPhase, + nb::arg("fraction"), nb::arg("coeffMarch"), + nb::arg("h"), nb::arg("k"), nb::arg("l")) + .def("SetPhasePar", &TextureMarchDollase::SetPhasePar, + nb::arg("i"), nb::arg("fraction"), nb::arg("coeffMarch"), + nb::arg("h"), nb::arg("k"), nb::arg("l")) + .def("DeletePhase", &TextureMarchDollase::DeletePhase, nb::arg("i")) + .def("GetNbPhase", &TextureMarchDollase::GetNbPhase) + .def("GetFraction", &TextureMarchDollase::GetFraction, nb::arg("i")) + .def("GetMarchCoeff", &TextureMarchDollase::GetMarchCoeff, nb::arg("i")) + .def("GetPhaseH", &TextureMarchDollase::GetPhaseH, nb::arg("i")) + .def("GetPhaseK", &TextureMarchDollase::GetPhaseK, nb::arg("i")) + .def("GetPhaseL", &TextureMarchDollase::GetPhaseL, nb::arg("i")) ; } diff --git a/src/extensions/nb_refinableobj.cpp b/src/extensions/nb_refinableobj.cpp index 4bdb782..322678d 100644 --- a/src/extensions/nb_refinableobj.cpp +++ b/src/extensions/nb_refinableobj.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include "helpers_nb.hpp" @@ -154,6 +155,8 @@ RefinableObj& extractRefinableObjArg(nb::object obj) try { return static_cast(nb::cast(obj)); } catch (...) {} try { return static_cast(nb::cast(obj)); } catch (...) {} try { return static_cast(nb::cast(obj)); } catch (...) {} + // Non-virtual, but RefinableObj is its second base (non-zero offset). + try { return static_cast(nb::cast(obj)); } catch (...) {} // Nothing matched -- let this raise a clear nanobind TypeError rather // than silently returning a dangling reference. return nb::cast(obj); diff --git a/src/extensions/nb_reflectionprofile.cpp b/src/extensions/nb_reflectionprofile.cpp index 14e442f..1252fd7 100644 --- a/src/extensions/nb_reflectionprofile.cpp +++ b/src/extensions/nb_reflectionprofile.cpp @@ -10,6 +10,7 @@ #undef B0 #include +#include #include #include "helpers_nb.hpp" @@ -133,6 +134,31 @@ void wrap_reflectionprofile(nb::module_& m) "Set the TCH pseudo-Voigt U, V, W, X, Y, Z, P and LGmix parameters.") ; + nb::class_( + m, "ReflectionProfileDoubleExponentialPseudoVoigt", + nb::is_weak_referenceable()) + .def(nb::init<>()) + .def(nb::init(), nb::arg("cell"), nb::keep_alive<1, 2>()) + .def(nb::init(), + nb::arg("old"), nb::keep_alive<1, 2>()) + .def("CreateCopy", + &ReflectionProfileDoubleExponentialPseudoVoigt::CreateCopy, + nb::rv_policy::take_ownership, nb::keep_alive<0, 1>(), + "Return an independent copy of this profile.") + .def("SetProfilePar", + &ReflectionProfileDoubleExponentialPseudoVoigt::SetProfilePar, + nb::arg("instrumentAlpha0"), nb::arg("instrumentAlpha1"), + nb::arg("instrumentBeta0"), nb::arg("instrumentBeta1"), + nb::arg("gaussianSigma0"), nb::arg("gaussianSigma1"), + nb::arg("gaussianSigma2"), nb::arg("lorentzianGamma0"), + nb::arg("lorentzianGamma1"), nb::arg("lorentzianGamma2"), + "Set the double-exponential pseudo-Voigt TOF profile parameters.") + .def("SetUnitCell", + &ReflectionProfileDoubleExponentialPseudoVoigt::SetUnitCell, + nb::arg("cell"), nb::keep_alive<1, 2>(), + "Set the unit cell used to compute d_hkl for reflections.") + ; + nb::class_( m, "ReflectionProfilePseudoVoigtAnisotropic") .def(nb::init<>()) diff --git a/src/pyobjcryst/powderpattern.py b/src/pyobjcryst/powderpattern.py index a044024..0b74b26 100644 --- a/src/pyobjcryst/powderpattern.py +++ b/src/pyobjcryst/powderpattern.py @@ -30,6 +30,7 @@ "PowderPatternDiffraction", "ReflectionProfileType", "SpaceGroupExplorer", + "TextureMarchDollase", ] from pyobjcryst._pyobjcryst import ( @@ -45,6 +46,7 @@ PowderPatternDiffraction, ReflectionProfileType, SpaceGroupExplorer, + TextureMarchDollase, ) from pyobjcryst.general import ObjCrystException from pyobjcryst.refinableobj import refpartype_scattdata_background diff --git a/src/pyobjcryst/reflectionprofile.py b/src/pyobjcryst/reflectionprofile.py index c3d9d3b..cf90ebd 100644 --- a/src/pyobjcryst/reflectionprofile.py +++ b/src/pyobjcryst/reflectionprofile.py @@ -43,6 +43,11 @@ aniso.SetProfilePar(fwhmCagliotiW=1e-4) ppd3.SetProfile(aniso) +For time-of-flight neutron data, use the double-exponential pseudo-Voigt +profile ``ReflectionProfileDoubleExponentialPseudoVoigt``. It needs the +unit cell to compute ``d_hkl``, which ``PowderPatternDiffraction.SetProfile`` +sets from the component's crystal. + Example ------- @@ -76,11 +81,13 @@ "ReflectionProfilePseudoVoigt", "ReflectionProfilePseudoVoigtTCH", "ReflectionProfilePseudoVoigtAnisotropic", + "ReflectionProfileDoubleExponentialPseudoVoigt", "ReflectionProfileType", ] from pyobjcryst._pyobjcryst import ( ReflectionProfile, + ReflectionProfileDoubleExponentialPseudoVoigt, ReflectionProfilePseudoVoigt, ReflectionProfilePseudoVoigtAnisotropic, ReflectionProfilePseudoVoigtTCH, diff --git a/tests/test_lsq.py b/tests/test_lsq.py index 0618d83..99730d2 100644 --- a/tests/test_lsq.py +++ b/tests/test_lsq.py @@ -93,18 +93,9 @@ def test_lsq_refine(self): self.c.RandomizeConfiguration() lsq.Refine(10, False, True) - def test_lsq_compiled_refined_obj_erase_all_param_set(self): - """Check parameter sets can be erased on compiled refined objects.""" - lsq = LSQ() - lsq.SetRefinedObj(self.d, 0, True, True) - lsq.PrepareRefParList() - refobj = lsq.GetCompiledRefinedObj() - save = refobj.CreateParamSet("save") - refobj.EraseAllParamSet() - self.assertRaises(ObjCrystException, refobj.SaveParamSet, save) - def test_lsq_get_variance_covariance_map(self): - """GetVarianceCovarianceMap returns a dict keyed by (name, name) tuples.""" + """GetVarianceCovarianceMap returns a dict keyed by (name, name) + tuples.""" lsq = LSQ() lsq.SetRefinedObj(self.d, 0, True, True) lsq.PrepareRefParList() @@ -137,7 +128,7 @@ def test_lsq_get_rw_history(self): self.assertGreaterEqual(v, 0.0) def test_lsq_refine_minrwpvar(self): - """minRwpVar argument is accepted by Refine without error.""" + """MinRwpVar argument is accepted by Refine without error.""" lsq = LSQ() lsq.SetRefinedObj(self.d, 0, True, True) lsq.PrepareRefParList() @@ -155,4 +146,3 @@ def test_lsq_prepare_ref_par_list_verbose(self): if __name__ == "__main__": unittest.main() - diff --git a/tests/test_powderpattern.py b/tests/test_powderpattern.py index 507734f..84fd164 100644 --- a/tests/test_powderpattern.py +++ b/tests/test_powderpattern.py @@ -23,7 +23,11 @@ from pyobjcryst import ObjCrystException from pyobjcryst.indexing import CrystalCentering, CrystalSystem, quick_index -from pyobjcryst.powderpattern import PowderPattern, SpaceGroupExplorer +from pyobjcryst.powderpattern import ( + PowderPattern, + SpaceGroupExplorer, + TextureMarchDollase, +) from pyobjcryst.radiation import RadiationType, WavelengthType from pyobjcryst.reflectionprofile import ReflectionProfileType @@ -365,16 +369,18 @@ def test_update_nbrefl(self): class _PowderPatternComponentFixtureMixin: - """Shared setUp for the three test classes below: a PowderPattern with - one real PowderPatternBackground and one real PowderPatternDiffraction - attached, mirroring TestPowderPattern.test_GetClassName's setup. Both - PowderPatternComponent and ScatteringData inherit RefinableObj virtually - in C++ (and PowderPatternDiffraction inherits PowderPatternComponent - virtually *again*, on top of ScatteringData non-virtually) -- see - nanobind_migration_notes.md's "Follow-up" section -- so neither class is - directly instantiable from Python; these fixtures exercise the - forwarding fix through real, concrete PowderPatternComponent-family - objects instead. + """Shared setUp for the three test classes below: a PowderPattern + with one real PowderPatternBackground and one real + PowderPatternDiffraction attached, mirroring + TestPowderPattern.test_GetClassName's setup. + + Both PowderPatternComponent and ScatteringData inherit RefinableObj + virtually in C++ (and PowderPatternDiffraction inherits + PowderPatternComponent virtually *again*, on top of ScatteringData + non-virtually) -- see nanobind_migration_notes.md's "Follow-up" + section -- so neither class is directly instantiable from Python; + these fixtures exercise the forwarding fix through real, concrete + PowderPatternComponent-family objects instead. """ @pytest.fixture(autouse=True) @@ -392,7 +398,9 @@ def setUp(self): self.diff = self.pp.AddPowderPatternDiffraction(c) -class TestPowderPatternComponent(_PowderPatternComponentFixtureMixin, unittest.TestCase): +class TestPowderPatternComponent( + _PowderPatternComponentFixtureMixin, unittest.TestCase +): def test_GetParentPowderPattern(self): self.assertIs(self.pp, self.bg.GetParentPowderPattern()) @@ -400,9 +408,10 @@ def test_GetParentPowderPattern(self): def test_refinableobj_forwarding(self): """PowderPatternComponent inherits RefinableObj virtually; only - GetName/SetName/GetClassName/FixAllPar/UnFixAllPar were reachable - before the forwarding fix (see test_GetClassName above, which - covers GetClassName specifically) -- exercise the rest here.""" + GetName/SetName/GetClassName/FixAllPar/UnFixAllPar were + reachable before the forwarding fix (see test_GetClassName + above, which covers GetClassName specifically) -- exercise the + rest here.""" self.bg.SetName("mybg") self.assertEqual("mybg", self.bg.GetName()) self.bg.FixAllPar() @@ -418,7 +427,9 @@ def test_refinableobj_forwarding(self): # ---------------------------------------------------------------------------- -class TestPowderPatternBackground(_PowderPatternComponentFixtureMixin, unittest.TestCase): +class TestPowderPatternBackground( + _PowderPatternComponentFixtureMixin, unittest.TestCase +): # def test___init__(self): assert False # def test_FixParametersBeyondMaxresolution(self): assert False # def test_GetPowderPatternCalc(self): assert False @@ -440,7 +451,9 @@ def test_refinableobj_forwarding(self): # ---------------------------------------------------------------------------- -class TestPowderPatternDiffraction(_PowderPatternComponentFixtureMixin, unittest.TestCase): +class TestPowderPatternDiffraction( + _PowderPatternComponentFixtureMixin, unittest.TestCase +): # def test___init__(self): assert False # def test_ExtractLeBail(self): assert False # def test_GetExtractionMode(self): assert False @@ -452,10 +465,10 @@ class TestPowderPatternDiffraction(_PowderPatternComponentFixtureMixin, unittest # def test_SetReflectionProfilePar(self): assert False def test_refinableobj_forwarding(self): - """PowderPatternDiffraction reaches RefinableObj via ScatteringData - (non-virtual, already a correct nanobind base) -- confirm the - forwarding applied to ScatteringData actually cascades here rather - than assuming it does.""" + """PowderPatternDiffraction reaches RefinableObj via + ScatteringData (non-virtual, already a correct nanobind base) -- + confirm the forwarding applied to ScatteringData actually + cascades here rather than assuming it does.""" self.diff.FixAllPar() self.diff.UnFixAllPar() self.assertEqual(0.0, self.diff.GetLogLikelihood()) @@ -497,5 +510,201 @@ def test_X2XCorrPhase_flat_detector_displacement(self): # ---------------------------------------------------------------------------- + +def build_textured_pattern(crystal, phases): + """Return an X-ray powder pattern of crystal with March-Dollase + phases.""" + pattern = PowderPattern() + pattern.SetWavelength(1.54) + x = np.deg2rad(np.linspace(5, 40, 2001)) + pattern.SetPowderPatternX(x) + pattern.SetPowderPatternObs(np.ones_like(x)) + diffraction = pattern.AddPowderPatternDiffraction(crystal) + texture = diffraction.GetTextureMarchDollase() + for phase in phases: + texture.AddPhase(*phase) + return pattern + + +@pytest.mark.parametrize( + "phases, expected_par_names", + [ + # Test that AddPhase stores a texture phase and its refinable + # parameters + # C1: One phase with a fraction, March coefficient and hkl is + # added. + # Expected: The getters return the phase's values, and its five + # parameters are named with the suffix _0. + ( + [(1.0, 0.5, 0, 0, 1)], + ["Fraction_0", "MarchCoeff_0", "H_0", "K_0", "L_0"], + ), + # C2: Two phases are added. + # Expected: The getters return each phase's values, and the second + # phase's parameters are named with the suffix _1. + ( + [(0.7, 0.5, 0, 0, 1), (0.3, 2.0, 1, 1, 0)], + [ + "Fraction_0", + "MarchCoeff_0", + "H_0", + "K_0", + "L_0", + "Fraction_1", + "MarchCoeff_1", + "H_1", + "K_1", + "L_1", + ], + ), + ], +) +def test_add_phase(loadcifdata, phases, expected_par_names): + crystal = loadcifdata("paracetamol.cif") + pattern = build_textured_pattern(crystal, phases) + texture = pattern.GetPowderPatternComponent(0).GetTextureMarchDollase() + actual_type = type(texture) + expected_type = TextureMarchDollase + actual_phases = [ + ( + texture.GetFraction(i), + texture.GetMarchCoeff(i), + texture.GetPhaseH(i), + texture.GetPhaseK(i), + texture.GetPhaseL(i), + ) + for i in range(texture.GetNbPhase()) + ] + expected_phases = phases + actual_par_names = [ + texture.GetPar(i).GetName() for i in range(texture.GetNbPar()) + ] + assert actual_type == expected_type + assert actual_phases == [pytest.approx(e) for e in expected_phases] + assert actual_par_names == expected_par_names + + +@pytest.mark.parametrize( + "initial_phases, method_name, method_args, expected_phases", + [ + # Test that changing the texture of an existing pattern + # recalculates it + # C1: A phase is added to an untextured pattern. + # Expected: The pattern equals one built with that phase. + ([], "AddPhase", (1.0, 0.5, 0, 0, 1), [(1.0, 0.5, 0, 0, 1)]), + # C2: SetPhasePar changes the March coefficient of a phase. + # Expected: The pattern equals one built with the new March + # coefficient. + ( + [(1.0, 0.5, 0, 0, 1)], + "SetPhasePar", + (0, 1.0, 2.0, 0, 0, 1), + [(1.0, 2.0, 0, 0, 1)], + ), + # C3: DeletePhase removes the second of two phases. + # Expected: The pattern equals one built with only the first phase. + ( + [(0.7, 0.5, 0, 0, 1), (0.3, 2.0, 1, 1, 0)], + "DeletePhase", + (1,), + [(0.7, 0.5, 0, 0, 1)], + ), + ], +) +def test_get_powder_pattern_calc_after_texture_change( + loadcifdata, initial_phases, method_name, method_args, expected_phases +): + crystal = loadcifdata("paracetamol.cif") + pattern = build_textured_pattern(crystal, initial_phases) + pattern.GetPowderPatternCalc() + texture = pattern.GetPowderPatternComponent(0).GetTextureMarchDollase() + getattr(texture, method_name)(*method_args) + actual_calc = pattern.GetPowderPatternCalc() + expected_calc = build_textured_pattern( + crystal, expected_phases + ).GetPowderPatternCalc() + assert np.allclose(actual_calc, expected_calc) + + +@pytest.mark.parametrize( + "method_name, method_args", + [ + # Test texture methods with the index of a phase that does not exist + # C1: SetPhasePar is called with index 1 on a single-phase texture. + # Expected: An ObjCrystException is raised. + ("SetPhasePar", (1, 1.0, 2.0, 0, 0, 1)), + # C2: DeletePhase is called with index 1 on a single-phase texture. + # Expected: An ObjCrystException is raised. + ("DeletePhase", (1,)), + ], +) +def test_texture_phase_index_out_of_range( + loadcifdata, method_name, method_args +): + crystal = loadcifdata("paracetamol.cif") + pattern = build_textured_pattern(crystal, [(1.0, 0.5, 0, 0, 1)]) + texture = pattern.GetPowderPatternComponent(0).GetTextureMarchDollase() + with pytest.raises(ObjCrystException): + getattr(texture, method_name)(*method_args) + + +@pytest.mark.parametrize( + "phases, deleted_index, expected_phases, expected_par_names", + [ + # Test that DeletePhase removes a phase and its refinable parameters + # C1: The only phase is deleted. + # Expected: No phase and no parameter remain. + ([(1.0, 0.5, 0, 0, 1)], 0, [], []), + # C2: The first of two phases is deleted. + # Expected: The second phase remains, and its parameters are + # renamed with the suffix _0. + ( + [(0.7, 0.5, 0, 0, 1), (0.3, 2.0, 1, 1, 0)], + 0, + [(0.3, 2.0, 1, 1, 0)], + ["Fraction_0", "MarchCoeff_0", "H_0", "K_0", "L_0"], + ), + ], +) +def test_delete_phase( + loadcifdata, phases, deleted_index, expected_phases, expected_par_names +): + crystal = loadcifdata("paracetamol.cif") + pattern = build_textured_pattern(crystal, phases) + texture = pattern.GetPowderPatternComponent(0).GetTextureMarchDollase() + texture.DeletePhase(deleted_index) + actual_phases = [ + ( + texture.GetFraction(i), + texture.GetMarchCoeff(i), + texture.GetPhaseH(i), + texture.GetPhaseK(i), + texture.GetPhaseL(i), + ) + for i in range(texture.GetNbPhase()) + ] + actual_par_names = [ + texture.GetPar(i).GetName() for i in range(texture.GetNbPar()) + ] + assert actual_phases == [pytest.approx(e) for e in expected_phases] + assert actual_par_names == expected_par_names + + +def test_get_powder_pattern_calc_after_march_coeff_par_change(loadcifdata): + # Test refining the March coefficient through its RefinablePar + # Expected: The pattern equals one built with the new March + # coefficient. + crystal = loadcifdata("paracetamol.cif") + pattern = build_textured_pattern(crystal, [(1.0, 0.5, 0, 0, 1)]) + pattern.GetPowderPatternCalc() + texture = pattern.GetPowderPatternComponent(0).GetTextureMarchDollase() + texture.GetPar("MarchCoeff_0").SetValue(2.0) + actual_calc = pattern.GetPowderPatternCalc() + expected_calc = build_textured_pattern( + crystal, [(1.0, 2.0, 0, 0, 1)] + ).GetPowderPatternCalc() + assert np.allclose(actual_calc, expected_calc) + + if __name__ == "__main__": unittest.main() diff --git a/tests/test_refinableobj.py b/tests/test_refinableobj.py index e5b71e5..087b3fe 100644 --- a/tests/test_refinableobj.py +++ b/tests/test_refinableobj.py @@ -425,7 +425,9 @@ def test_xml(self): def test_xml_output_input_reachable(self): """XMLOutput/XMLInput were dropped from the nanobind port's - method-forwarding audit alongside the LSQ methods below. Note + method-forwarding audit alongside the LSQ methods below. + + Note RefinableObj's own C++ base implementation is a documented no-op ("Does nothing! Should be purely virtual...", see RefinableObj/IO.cpp) -- real content only appears once a concrete @@ -434,7 +436,8 @@ def test_xml_output_input_reachable(self): file-like object, exercised on concrete subclasses that do override XMLOutput/XMLInput with real behaviour, live in test_scatterer.py (Atom), test_powderpattern.py (PowderPatternBackground) and - test_single_crystal_data.py (DiffractionDataSingleCrystal).""" + test_single_crystal_data.py (DiffractionDataSingleCrystal). + """ import io buf = io.StringIO() diff --git a/tests/test_reflectionprofile.py b/tests/test_reflectionprofile.py index f4cfc48..a3b1eeb 100644 --- a/tests/test_reflectionprofile.py +++ b/tests/test_reflectionprofile.py @@ -1,13 +1,17 @@ """Unit tests for pyobjcryst.reflectionprofile bindings.""" +import gc import unittest +import weakref import numpy as np import pytest from pyobjcryst.powderpattern import PowderPattern +from pyobjcryst.radiation import RadiationType, WavelengthType from pyobjcryst.refinableobj import RefinableObj from pyobjcryst.reflectionprofile import ( + ReflectionProfileDoubleExponentialPseudoVoigt, ReflectionProfilePseudoVoigt, ReflectionProfilePseudoVoigtAnisotropic, ReflectionProfilePseudoVoigtTCH, @@ -126,8 +130,8 @@ def test_is_anisotropic(self): self.assertFalse(self.profile.IsAnisotropic()) def test_concrete_pseudo_voigt_profiles(self): - """ReflectionProfilePseudoVoigt and Anisotropic can be constructed - and installed on a PowderPatternDiffraction.""" + """ReflectionProfilePseudoVoigt and Anisotropic can be + constructed and installed on a PowderPatternDiffraction.""" iso = ReflectionProfilePseudoVoigt() self.assertFalse(iso.IsAnisotropic()) self.assertTrue( @@ -148,7 +152,9 @@ def test_concrete_pseudo_voigt_profiles(self): aniso.SetProfilePar(fwhmCagliotiW=1e-4) self.ppd.SetProfile(aniso) installed_aniso = self.ppd.GetProfile() - self.assertIsInstance(installed_aniso, ReflectionProfilePseudoVoigtAnisotropic) + self.assertIsInstance( + installed_aniso, ReflectionProfilePseudoVoigtAnisotropic + ) def test_tch_profile_width_and_limits(self): """TCH derives one common FWHM and reaches both pure limits.""" @@ -230,5 +236,140 @@ def test_anisotropic_set_profile_par(self): self.assertAlmostEqual(aniso.GetPar("G_KK").GetValue(), 2e-5) +def build_tof_pattern(crystal, profile): + """Return a TOF neutron powder pattern of crystal using profile.""" + pattern = PowderPattern() + pattern.SetRadiationType(RadiationType.RAD_NEUTRON) + pattern.GetRadiation().SetWavelengthType(WavelengthType.WAVELENGTH_TOF) + # TOF decreases with sin(theta)/lambda, so x is in descending order. + x = np.linspace(120000, 20000, 2001) + pattern.SetPowderPatternX(x) + pattern.SetPowderPatternObs(np.ones_like(x)) + diffraction = pattern.AddPowderPatternDiffraction(crystal) + diffraction.SetProfile(profile) + return pattern + + +def test_set_profile_par_tof(): + # Test that SetProfilePar assigns the ten TOF profile parameters + # Expected: Each refinable parameter has the value passed for it. + profile = ReflectionProfileDoubleExponentialPseudoVoigt() + profile.SetProfilePar( + instrumentAlpha0=0.01, + instrumentAlpha1=0.1, + instrumentBeta0=0.02, + instrumentBeta1=0.004, + gaussianSigma0=1.0, + gaussianSigma1=7.0, + gaussianSigma2=0.5, + lorentzianGamma0=0.1, + lorentzianGamma1=0.2, + lorentzianGamma2=0.4, + ) + actual_values = { + profile.GetPar(i).GetName(): profile.GetPar(i).GetValue() + for i in range(profile.GetNbPar()) + } + expected_values = { + "Alpha0": 0.01, + "Alpha1": 0.1, + "Beta0": 0.02, + "Beta1": 0.004, + "GaussianSigma0": 1.0, + "GaussianSigma1": 7.0, + "GaussianSigma2": 0.5, + "LorentzianGamma0": 0.1, + "LorentzianGamma1": 0.2, + "LorentzianGamma2": 0.4, + } + assert actual_values == pytest.approx(expected_values) + + +def test_set_profile_tof_without_unit_cell(loadcifdata): + # Test installing a TOF profile constructed without a unit cell + # Expected: SetProfile gives the profile the crystal's unit cell, so + # the strongest peak is the Si (220) reflection at + # TOF = DIFC * d + DIFA * d^2. + crystal = loadcifdata("Si.cif") + pattern = build_tof_pattern( + crystal, ReflectionProfileDoubleExponentialPseudoVoigt() + ) + calc = pattern.GetPowderPatternCalc() + x = pattern.GetPowderPatternX() + actual_peak_tof = x[np.argmax(calc)] + d_220 = crystal.a / np.sqrt(8) + expected_peak_tof = ( + pattern.GetPar("DIFC").GetValue() * d_220 + + pattern.GetPar("DIFA").GetValue() * d_220**2 + ) + tof_step = abs(x[1] - x[0]) + assert actual_peak_tof == pytest.approx(expected_peak_tof, abs=tof_step) + + +@pytest.mark.parametrize( + "parameter_name, new_value", + [ + # Test that changing a parameter of the installed TOF profile + # recalculates the pattern + # C1: The rising-edge exponential coefficient Alpha1 is changed. + # Expected: The pattern equals one calculated with Alpha1 set + # before the profile was installed. + ("Alpha1", 0.2), + # C2: The falling-edge exponential coefficient Beta0 is changed. + # Expected: The pattern equals one calculated with Beta0 set + # before the profile was installed. + ("Beta0", 0.05), + # C3: The Gaussian width coefficient GaussianSigma1 is changed. + # Expected: The pattern equals one calculated with GaussianSigma1 + # set before the profile was installed. + ("GaussianSigma1", 50.0), + # C4: The Lorentzian width coefficient LorentzianGamma2 is changed. + # Expected: The pattern equals one calculated with + # LorentzianGamma2 set before the profile was installed. + ("LorentzianGamma2", 2.0), + ], +) +def test_get_powder_pattern_calc_after_tof_profile_change( + loadcifdata, parameter_name, new_value +): + crystal = loadcifdata("Si.cif") + pattern = build_tof_pattern( + crystal, ReflectionProfileDoubleExponentialPseudoVoigt() + ) + pattern.GetPowderPatternCalc() + profile = pattern.GetPowderPatternComponent(0).GetProfile() + profile.GetPar(parameter_name).SetValue(new_value) + actual_calc = pattern.GetPowderPatternCalc() + expected_profile = ReflectionProfileDoubleExponentialPseudoVoigt(crystal) + expected_profile.GetPar(parameter_name).SetValue(new_value) + expected_calc = build_tof_pattern( + crystal, expected_profile + ).GetPowderPatternCalc() + assert np.allclose(actual_calc, expected_calc) + + +@pytest.mark.parametrize("copy_method", ["CreateCopy", "copy_constructor"]) +def test_tof_profile_copy_keeps_original_alive(loadcifdata, copy_method): + # Test that a copy of a TOF profile keeps its source and unit cell alive + # Expected: The original profile and crystal survive until the copy is + # deleted, since both profiles reference the same unit cell. + crystal = loadcifdata("Si.cif") + original = ReflectionProfileDoubleExponentialPseudoVoigt(crystal) + if copy_method == "CreateCopy": + copied = original.CreateCopy() + else: + copied = ReflectionProfileDoubleExponentialPseudoVoigt(original) + original_ref = weakref.ref(original) + crystal_ref = weakref.ref(crystal) + del original, crystal + gc.collect() + assert original_ref() is not None + assert crystal_ref() is not None + del copied + gc.collect() + assert original_ref() is None + assert crystal_ref() is None + + if __name__ == "__main__": unittest.main() diff --git a/tests/test_scatterer.py b/tests/test_scatterer.py index c68c89c..d2abece 100644 --- a/tests/test_scatterer.py +++ b/tests/test_scatterer.py @@ -1,10 +1,11 @@ #!/usr/bin/env python """Functional tests for the RefinableObj/ScatteringPower method surface restored on Scatterer/ScatteringPowerAtom (see -nanobind_migration_notes.md's "Follow-up: restoring RefinableObj's method -surface on virtually-inherited classes"; presence-only coverage across every -affected class lives in test_virtualbase_forwarding.py). These exercise the -actual behaviour, not just that the method exists. +nanobind_migration_notes.md's "Follow-up: restoring RefinableObj's +method surface on virtually-inherited classes"; presence-only coverage +across every affected class lives in test_virtualbase_forwarding.py). + +These exercise the actual behaviour, not just that the method exists. """ import io @@ -17,7 +18,8 @@ class TestScattererForwarding(unittest.TestCase): - """Atom reaches RefinableObj only through Scatterer's C++ virtual base.""" + """Atom reaches RefinableObj only through Scatterer's C++ virtual + base.""" def setUp(self): self.sp, self.atom = makeScatterer() @@ -62,9 +64,10 @@ def test_update_display(self): self.atom.UpdateDisplay() def test_get_clock_scatt_comp_list(self): - """Forwards a *protected* Scatterer member (mClockScattCompList); - unrelated to the virtual-inheritance issue (it's a plain omission) - but fixed in the same pass -- see PyScatterer in nb_scatterer.cpp.""" + """Forwards a *protected* Scatterer member + (mClockScattCompList); unrelated to the virtual-inheritance + issue (it's a plain omission) but fixed in the same pass -- see + PyScatterer in nb_scatterer.cpp.""" clock = self.atom.GetClockScattCompList() self.assertIsInstance(clock, RefinableObjClock) @@ -86,9 +89,10 @@ def test_create_param_set_roundtrip(self): class TestScatteringPowerAtomForwarding(unittest.TestCase): """ScatteringPowerAtom inherits ScatteringPower virtually, and - ScatteringPower inherits RefinableObj virtually too -- a second hop, so - both bind_refinableobj_forwarding() and - bind_scatteringpower_forwarding() had to be applied directly to it.""" + ScatteringPower inherits RefinableObj virtually too -- a second hop, + so both bind_refinableobj_forwarding() and + bind_scatteringpower_forwarding() had to be applied directly to + it.""" def setUp(self): self.sp, self.atom = makeScatterer() @@ -105,8 +109,9 @@ def test_log_likelihood(self): self.assertEqual(0.0, self.sp.GetLogLikelihood()) def test_scatteringpower_own_methods_forwarded(self): - """These are ScatteringPower's *own* methods (not RefinableObj's), - forwarded via bind_scatteringpower_forwarding().""" + """These are ScatteringPower's *own* methods (not + RefinableObj's), forwarded via + bind_scatteringpower_forwarding().""" self.assertGreaterEqual(self.sp.GetDynPopCorrIndex(), 0) self.assertIsNotNone(self.sp.GetLastChangeClock()) r, g, b = self.sp.GetColourRGB() @@ -114,8 +119,9 @@ def test_scatteringpower_own_methods_forwarded(self): def test_lsq_methods_reachable_and_raise_when_unset(self): """Ported forward alongside the forwarding fix: GetLSQCalc/Obs/ - Weight/Deriv exist on RefinableObj and are therefore forwarded here - too, even though nothing has configured a custom LSQ target.""" + Weight/Deriv exist on RefinableObj and are therefore forwarded + here too, even though nothing has configured a custom LSQ + target.""" with self.assertRaises(ObjCrystException): self.sp.GetLSQCalc(0) diff --git a/tests/test_virtualbase_forwarding.py b/tests/test_virtualbase_forwarding.py index c497f29..be32f2b 100644 --- a/tests/test_virtualbase_forwarding.py +++ b/tests/test_virtualbase_forwarding.py @@ -1,5 +1,6 @@ #!/usr/bin/env python -"""Coverage guard for the RefinableObj/ScatteringPower method-forwarding fix. +"""Coverage guard for the RefinableObj/ScatteringPower method-forwarding +fix. Scatterer, ScatteringPower, ScatteringData, and PowderPatternComponent all inherit RefinableObj *virtually* in C++, so nanobind cannot declare @@ -46,8 +47,8 @@ def _public_methods(cls): class TestRefinableObjSurfaceForwarding(unittest.TestCase): - """Every RefinableObj method must be reachable on the four classes that - only reach RefinableObj through a C++ virtual base.""" + """Every RefinableObj method must be reachable on the four classes + that only reach RefinableObj through a C++ virtual base.""" def _assert_full_coverage(self, cls): missing = [m for m in REFINABLEOBJ_METHODS if not hasattr(cls, m)] @@ -108,9 +109,10 @@ def test_globalscatteringpower(self): class TestScatteringPowerSurfaceForwarding(unittest.TestCase): - """ScatteringPowerAtom and GlobalScatteringPower need ScatteringPower's - own method surface (GetScatteringFactor, GetBij/SetBij, GetColour, ...) - forwarded too, on top of RefinableObj's.""" + """ScatteringPowerAtom and GlobalScatteringPower need + ScatteringPower's own method surface (GetScatteringFactor, + GetBij/SetBij, GetColour, ...) forwarded too, on top of + RefinableObj's.""" def _assert_full_coverage(self, cls): missing = [m for m in SCATTERINGPOWER_METHODS if not hasattr(cls, m)]