DiffractionDataSingleCrystal: add SetHklIobs, SetIobs, SetSigma, GetSigma, GetChi2, FitScaleFactorForRw and FitScaleFactorForR. Add a unittest.

Fixes #42

Author: vincefn

src/extensions/diffractiondatasinglecrystal_ext.cpp

@@ -13,13 +13,17 @@
 ******************************************************************************
 *
 * boost::python bindings to ObjCryst::DiffractionDataSingleCrystal.
-*
+* Changes from ObjCryst::DiffractionDataSingleCrystal
+* - SetHklIobs takes float(64) H, K, L arrays rather than long integers - easier
+* because passing numpy int array seesm complicated, and more practical anyway
+* since GetH() GetK() GetL() functions from ScatteringData natively use floats.
 *****************************************************************************/
 
 #include <boost/python/class.hpp>
 #include <boost/python/def.hpp>
 #include <boost/python/copy_const_reference.hpp>
 #include <boost/python/manage_new_object.hpp>
+#include <boost/python/stl_iterator.hpp>
 #undef B0
 
 #include <string>
@@ -67,6 +71,71 @@ DiffractionDataSingleCrystal* _CreateSingleCrystalDataFromCIF(bp::object input,
     return d;
 }
 
+void setdiffractiondatasinglecrystal_iobs(DiffractionDataSingleCrystal& diff, bp::object iobs)
+{
+    CrystVector_REAL iiobs;
+    assignCrystVector(iiobs, iobs);
+    if(iiobs.size() != diff.GetIobs().size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetIobs(): "
+                                        "number of elements does not match the previous Iobs list. "
+                                        "Use SetHklIobs if you want to change the number of reflections.");
+    MuteObjCrystUserInfo muzzle;
+    diff.SetIobs(iiobs);
+}
+
+void setdiffractiondatasinglecrystal_sigma(DiffractionDataSingleCrystal& diff, bp::object sigma)
+{
+    CrystVector_REAL ssigma;
+    assignCrystVector(ssigma, sigma);
+    if(ssigma.size() != diff.GetIobs().size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetSigma(): "
+                                        "number of elements does not match the Iobs list. "
+                                        "Use SetHklIobs if you want to change the number of reflections.");
+    MuteObjCrystUserInfo muzzle;
+    diff.SetSigma(ssigma);
+}
+
+// TODO: For SetHklIobs we should pass directly an integer array but that seems difficult-passed numpy arrays
+//       are always interpreted as doubles (?). It's more practical this way.
+void assignCrystVector(CrystVector<long>& cv, bp::object obj)
+{
+    bp::stl_input_iterator<double> begin(obj), end;
+    std::list<double> values(begin, end);
+    cv.resize(values.size());
+    std::list<double>::const_iterator vv = values.begin();
+    long* dst = cv.data();
+    for (; vv != values.end(); ++vv, ++dst) *dst = lround(*vv);
+}
+
+void setdiffractiondatasinglecrystal_hkliobs(DiffractionDataSingleCrystal& diff,
+                                             bp::object h,bp::object k, bp::object l,
+                                             bp::object iobs, bp::object sigma)
+{
+    CrystVector<long> hh;
+    assignCrystVector(hh, h);
+    CrystVector<long> kk;
+    assignCrystVector(kk, k);
+    CrystVector<long> ll;
+    assignCrystVector(ll, l);
+    CrystVector_REAL iiobs;
+    assignCrystVector(iiobs, iobs);
+    CrystVector_REAL ssigma;
+    assignCrystVector(ssigma, sigma);
+
+    if(hh.size() != kk.size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetHklIobs(): h and k array sizes differ");
+    if(hh.size() != ll.size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetHklIobs(): h and l array sizes differ");
+    if(hh.size() != iiobs.size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetHklIobs(): h and iobs array sizes differ");
+    if(hh.size() != ssigma.size())
+      throw ObjCryst::ObjCrystException("DiffractionDataSingleCrystal::SetHklIobs(): h and sigma array sizes differ");
+
+    MuteObjCrystUserInfo muzzle;
+    diff.SetHklIobs(hh, kk, ll, iiobs, ssigma);
+}
+
+
 }   // namespace
 
 void wrap_diffractiondatasinglecrystal()
@@ -83,10 +152,24 @@ void wrap_diffractiondatasinglecrystal()
                 return_value_policy<copy_const_reference>())
         .def("GetIobs", &DiffractionDataSingleCrystal::GetIobs,
                 return_value_policy<copy_const_reference>())
-        // FIXME ... add SetIobs, SetSigma ....
+        .def("GetSigma", &DiffractionDataSingleCrystal::GetSigma,
+                return_value_policy<copy_const_reference>())
+        .def("SetIobs",
+                &setdiffractiondatasinglecrystal_iobs,
+                bp::arg("iobs"))
+        .def("SetSigma",
+                &setdiffractiondatasinglecrystal_sigma,
+                bp::arg("sigma"))
+        .def("SetHklIobs",
+                &setdiffractiondatasinglecrystal_hkliobs,
+                (bp::arg("h"),bp::arg("k"),bp::arg("l"), bp::arg("iobs"), bp::arg("sigma")))
         .def("SetIobsToIcalc", &DiffractionDataSingleCrystal::SetIobsToIcalc)
         .def("GetRw", &DiffractionDataSingleCrystal::GetRw)
         .def("GetR", &DiffractionDataSingleCrystal::GetR)
+        .def("GetChi2", &DiffractionDataSingleCrystal::GetChi2)
+        .def("FitScaleFactorForRw", &DiffractionDataSingleCrystal::FitScaleFactorForRw)
+        .def("FitScaleFactorForR", &DiffractionDataSingleCrystal::FitScaleFactorForR)
+        // TODO: These functions should print a limited number of reflections - problems otherwise
         .def("PrintObsData", &DiffractionDataSingleCrystal::PrintObsData)
         .def("PrintObsCalcData", &DiffractionDataSingleCrystal::PrintObsCalcData)
         .def("SetUseOnlyLowAngleData", &DiffractionDataSingleCrystal::SetUseOnlyLowAngleData)

src/pyobjcryst/tests/test_single_crystal_data.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+
+"""Tests for diffractiondatasinglecrystal module."""
+
+import unittest
+import gc
+import numpy as np
+
+from pyobjcryst.crystal import CreateCrystalFromCIF, Crystal
+from pyobjcryst.diffractiondatasinglecrystal import *
+from pyobjcryst.tests.pyobjcrysttestutils import loadcifdata, datafile
+
+
+class test_single_crystal_data(unittest.TestCase):
+
+    def test_create(self):
+        """Test creating a DiffractionDataSingleCrystal object"""
+        c = Crystal(3.52, 3.52, 3.52, "225")
+        d = DiffractionDataSingleCrystal(c)
+
+    def test_create_set_hkliobs(self):
+        """test SetHklIobs, SetIobs and SetSigma"""
+        c = Crystal(3.1, 3.2, 3.3, "Pmmm")
+        d = DiffractionDataSingleCrystal(c)
+        n0 = 5
+        nb = n0 ** 3
+        r = np.arange(1, nb + 1, dtype=np.float64)
+        h = r % n0
+        l = r // n0 ** 2
+        k = (r - l * n0 ** 2) // n0
+        iobs = np.random.uniform(0, 100, nb)
+        sigma = np.sqrt(iobs)
+
+        d.SetHklIobs(h, k, l, iobs, sigma)
+
+        # SetHklIobs sorts reflecions by sin(theta)/lambda, so do the same for comparison
+        s = np.sqrt(h ** 2 / 3.1 ** 2 + k ** 2 / 3.2 ** 2 + l ** 2 / 3.3 ** 2) / 2
+        idx = np.argsort(s)
+
+        iobs = np.take(iobs, idx)
+        sigma = np.take(sigma, idx)
+        h = np.take(h, idx)
+        k = np.take(k, idx)
+        l = np.take(l, idx)
+        self.assertTrue(np.all(iobs == d.GetIobs()))
+        self.assertTrue(np.all(sigma == d.GetSigma()))
+        self.assertTrue(np.all(h == d.GetH()))
+        self.assertTrue(np.all(k == d.GetK()))
+        self.assertTrue(np.all(l == d.GetL()))
+
+        # Set Iobs and sigma individually
+        iobs = np.random.uniform(0, 100, nb)
+        d.SetIobs(iobs)
+        self.assertTrue(np.all(iobs == d.GetIobs()))
+
+        sigma = np.random.uniform(0, 10, nb)
+        d.SetSigma(sigma)
+        self.assertTrue(np.all(sigma == d.GetSigma()))
+
+if __name__ == "__main__":
+    unittest.main()