feat: port March-Dollase texture (#119) and TOF profile (#120) to nanobind

.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) }}

.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:

.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:

.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:

.pre-commit-config.yaml

@@ -1,3 +1,5 @@
+# Vendored third-party sources (objcryst submodule, cctbx, newmat)
+exclude: "^extern/"
 default_language_version:
   python: python3
 ci:

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:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

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:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

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,

src/extensions/nb_powderpatterndiffraction.cpp

@@ -11,13 +11,46 @@
 #include <ObjCryst/ObjCryst/PowderPattern.h>
 #include <ObjCryst/ObjCryst/ScatteringData.h>
 #include <ObjCryst/ObjCryst/ReflectionProfile.h>
+#include <ObjCryst/ObjCryst/ScatteringCorr.h>
 #include <ObjCryst/RefinableObj/RefinableObj.h>
 
 #include "helpers_nb.hpp"
 
 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<ReflectionProfileDoubleExponentialPseudoVoigt*>(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<RefinableObj>& subobjs = pdiff.GetSubObjRegistry();
+    for (long i = 0; i < subobjs.GetNb(); ++i)
+    {
+        auto* texture = dynamic_cast<TextureMarchDollase*>(&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_<ReflectionProfileType>(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_<TextureMarchDollase>(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"))
         ;
 }

src/extensions/nb_refinableobj.cpp

@@ -28,6 +28,7 @@
 #include <ObjCryst/ObjCryst/PowderPattern.h>
 #include <ObjCryst/ObjCryst/Scatterer.h>
 #include <ObjCryst/ObjCryst/ScatteringPower.h>
+#include <ObjCryst/ObjCryst/ScatteringCorr.h>
 
 #include "helpers_nb.hpp"
 
@@ -154,6 +155,8 @@ RefinableObj& extractRefinableObjArg(nb::object obj)
     try { return static_cast<RefinableObj&>(nb::cast<Scatterer&>(obj)); } catch (...) {}
     try { return static_cast<RefinableObj&>(nb::cast<ScatteringPower&>(obj)); } catch (...) {}
     try { return static_cast<RefinableObj&>(nb::cast<ScatteringPowerAtom&>(obj)); } catch (...) {}
+    // Non-virtual, but RefinableObj is its second base (non-zero offset).
+    try { return static_cast<RefinableObj&>(nb::cast<TextureMarchDollase&>(obj)); } catch (...) {}
     // Nothing matched -- let this raise a clear nanobind TypeError rather
     // than silently returning a dangling reference.
     return nb::cast<RefinableObj&>(obj);

src/extensions/nb_reflectionprofile.cpp

@@ -10,6 +10,7 @@
 
 #undef B0
 #include <ObjCryst/ObjCryst/ReflectionProfile.h>
+#include <ObjCryst/ObjCryst/UnitCell.h>
 #include <ObjCryst/RefinableObj/RefinableObj.h>
 
 #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_<ReflectionProfileDoubleExponentialPseudoVoigt, ReflectionProfile>(
+        m, "ReflectionProfileDoubleExponentialPseudoVoigt",
+        nb::is_weak_referenceable())
+        .def(nb::init<>())
+        .def(nb::init<const UnitCell&>(), nb::arg("cell"), nb::keep_alive<1, 2>())
+        .def(nb::init<const ReflectionProfileDoubleExponentialPseudoVoigt&>(),
+             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_<ReflectionProfilePseudoVoigtAnisotropic, ReflectionProfile>(
         m, "ReflectionProfilePseudoVoigtAnisotropic")
         .def(nb::init<>())

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

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,

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()
-