From c981f715515883b128ea926fa3411b59be2712ce Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 15:36:16 +0100 Subject: [PATCH 01/26] Add a probe enum --- src/dxtbx/model/beam.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 9903341ec..beae21fca 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -24,6 +24,9 @@ namespace dxtbx { namespace model { using scitbx::vec3; + // probe type enumeration + enum probe { xray = 1, electron = 2, neutron = 3 }; + /** Base class for beam objects */ class BeamBase { public: @@ -385,6 +388,7 @@ namespace dxtbx { namespace model { double flux_; double transmission_; scitbx::af::shared > s0_at_scan_points_; + probe probe_; }; /** Print beam information */ From bae2822adddf66ffaa1c4b88926f7ca093d14585 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 15:42:24 +0100 Subject: [PATCH 02/26] get_probe accessor --- src/dxtbx/model/beam.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index beae21fca..4d49b01cf 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -47,6 +47,7 @@ namespace dxtbx { namespace model { virtual std::size_t get_num_scan_points() const = 0; virtual scitbx::af::shared > get_s0_at_scan_points() const = 0; virtual vec3 get_s0_at_scan_point(std::size_t index) const = 0; + virtual probe get_probe() const = 0; virtual void set_direction(vec3 direction) = 0; virtual void set_wavelength(double wavelength) = 0; @@ -292,6 +293,10 @@ namespace dxtbx { namespace model { return s0_at_scan_points_[index]; } + probe get_probe() const { + return probe_; + } + void reset_scan_points() { s0_at_scan_points_.clear(); } From 6eea2fd78a43eb0cdc87be2c0253fedaaf71455c Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 16:03:14 +0100 Subject: [PATCH 03/26] probe-->Probe, and export to Python. So now, can do from dxtbx.model.beam import Probe --- src/dxtbx/model/beam.h | 8 ++++---- src/dxtbx/model/beam.py | 4 ++-- src/dxtbx/model/boost_python/beam.cc | 6 ++++++ 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 4d49b01cf..0769f53ed 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -25,7 +25,7 @@ namespace dxtbx { namespace model { using scitbx::vec3; // probe type enumeration - enum probe { xray = 1, electron = 2, neutron = 3 }; + enum Probe { xray = 1, electron = 2, neutron = 3 }; /** Base class for beam objects */ class BeamBase { @@ -47,7 +47,7 @@ namespace dxtbx { namespace model { virtual std::size_t get_num_scan_points() const = 0; virtual scitbx::af::shared > get_s0_at_scan_points() const = 0; virtual vec3 get_s0_at_scan_point(std::size_t index) const = 0; - virtual probe get_probe() const = 0; + virtual Probe get_probe() const = 0; virtual void set_direction(vec3 direction) = 0; virtual void set_wavelength(double wavelength) = 0; @@ -293,7 +293,7 @@ namespace dxtbx { namespace model { return s0_at_scan_points_[index]; } - probe get_probe() const { + Probe get_probe() const { return probe_; } @@ -393,7 +393,7 @@ namespace dxtbx { namespace model { double flux_; double transmission_; scitbx::af::shared > s0_at_scan_points_; - probe probe_; + Probe probe_; }; /** Print beam information */ diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 870dd9c4b..f128e85a1 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -7,9 +7,9 @@ import libtbx.phil try: - from ..dxtbx_model_ext import Beam + from ..dxtbx_model_ext import Beam, Probe except ModuleNotFoundError: - from dxtbx_model_ext import Beam # type: ignore + from dxtbx_model_ext import Beam, Probe # type: ignore beam_phil_scope = libtbx.phil.parse( """ diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index bf709d526..a9b77751a 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -17,6 +17,7 @@ #include #include #include +#include namespace dxtbx { namespace model { namespace boost_python { namespace beam_detail { @@ -224,6 +225,11 @@ namespace dxtbx { namespace model { namespace boost_python { void export_beam() { using namespace beam_detail; + enum_("Probe") + .value("xray", xray) + .value("electron", electron) + .value("neutron", neutron); + class_("BeamBase", no_init) .def("get_sample_to_source_direction", &BeamBase::get_sample_to_source_direction) .def("set_direction", &BeamBase::set_direction) From eb22c37cc0852f9edd3105050fc4ee28d27c7772 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 16:08:40 +0100 Subject: [PATCH 04/26] Export get_probe Note the value is uninitialized! --- src/dxtbx/model/boost_python/beam.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index a9b77751a..f528d8902 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -260,6 +260,7 @@ namespace dxtbx { namespace model { namespace boost_python { .def("set_s0_at_scan_points", &Beam_set_s0_at_scan_points_from_list) .def("get_s0_at_scan_points", &BeamBase::get_s0_at_scan_points) .def("get_s0_at_scan_point", &BeamBase::get_s0_at_scan_point) + .def("get_probe", &BeamBase::get_probe) .def("reset_scan_points", &BeamBase::reset_scan_points) .def("rotate_around_origin", &rotate_around_origin, From 2b79fd8f6dbb72f351549b072421c66d548e9906 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 16:15:48 +0100 Subject: [PATCH 05/26] Ensure probe_ is initialized. Now we get: >>> from dxtbx.model.beam import Beam >>> b=Beam() >>> b.get_probe() dxtbx_model_ext.Probe.xray --- src/dxtbx/model/beam.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 0769f53ed..876d482cf 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -86,7 +86,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) {} + transmission_(1.0), + probe_(Probe::xray) {} /** * @param s0 The incident beam vector. From 6979e25aecf1a7e03ffd1884f443a4daabb670be Mon Sep 17 00:00:00 2001 From: David Waterman Date: Fri, 21 Jul 2023 16:34:19 +0100 Subject: [PATCH 06/26] Add set_probe and export. Now this works: >>> from dxtbx.model.beam import Beam >>> b=Beam() >>> b.get_probe() dxtbx_model_ext.Probe.xray >>> from dxtbx.model.beam import Probe >>> b.set_probe(Probe.neutron) >>> b.get_probe() dxtbx_model_ext.Probe.neutron --- src/dxtbx/model/beam.h | 5 +++++ src/dxtbx/model/boost_python/beam.cc | 1 + 2 files changed, 6 insertions(+) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 876d482cf..406b6ed52 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -63,6 +63,7 @@ namespace dxtbx { namespace model { virtual void set_transmission(double transmission) = 0; virtual void set_s0_at_scan_points( const scitbx::af::const_ref > &s0) = 0; + virtual void set_probe(Probe probe) = 0; virtual void reset_scan_points() = 0; virtual bool is_similar_to(const BeamBase &rhs, @@ -298,6 +299,10 @@ namespace dxtbx { namespace model { return probe_; } + void set_probe(Probe probe) { + probe_ = probe; + } + void reset_scan_points() { s0_at_scan_points_.clear(); } diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index f528d8902..370a22b28 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -261,6 +261,7 @@ namespace dxtbx { namespace model { namespace boost_python { .def("get_s0_at_scan_points", &BeamBase::get_s0_at_scan_points) .def("get_s0_at_scan_point", &BeamBase::get_s0_at_scan_point) .def("get_probe", &BeamBase::get_probe) + .def("set_probe", &BeamBase::set_probe) .def("reset_scan_points", &BeamBase::reset_scan_points) .def("rotate_around_origin", &rotate_around_origin, From e5c88399296552d07ebb184684705eb331b5af52 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 10:02:21 +0100 Subject: [PATCH 07/26] Ensure probe is initialised by the other constructors too --- src/dxtbx/model/beam.h | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 406b6ed52..ac161411f 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -99,7 +99,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(s0.length() > 0); wavelength_ = 1.0 / s0.length(); direction_ = -s0.normalize(); @@ -116,7 +117,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } @@ -132,7 +134,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(s0.length() > 0); wavelength_ = 1.0 / s0.length(); direction_ = -s0.normalize(); @@ -154,7 +157,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } @@ -183,7 +187,8 @@ namespace dxtbx { namespace model { polarization_normal_(polarization_normal), polarization_fraction_(polarization_fraction), flux_(flux), - transmission_(transmission) { + transmission_(transmission), + probe_(Probe::xray) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } From 48c096df385c2835984a52c8ec1c4484e86e1634 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 10:16:45 +0100 Subject: [PATCH 08/26] Incorporate probe into equality and similarity checks --- src/dxtbx/model/beam.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index ac161411f..55d0b603a 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -343,7 +343,8 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= eps && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= eps; + <= eps + && (probe_ == rhs.get_probe()); } bool is_similar_to(const BeamBase &rhs, @@ -380,7 +381,8 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= polarization_normal_tolerance && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= polarization_fraction_tolerance; + <= polarization_fraction_tolerance + && (probe_ == rhs.get_probe()); } bool operator!=(const BeamBase &rhs) const { From 899b166748c520281f85916423f1b9f10b64fae9 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 11:03:53 +0100 Subject: [PATCH 09/26] Add get_probe_name method so that the NeXus probe name can be returned for each of the defined probe values --- src/dxtbx/model/beam.h | 18 ++++++++++++++++++ src/dxtbx/model/boost_python/beam.cc | 1 + 2 files changed, 19 insertions(+) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 55d0b603a..f3f1ed6a2 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -48,6 +49,7 @@ namespace dxtbx { namespace model { virtual scitbx::af::shared > get_s0_at_scan_points() const = 0; virtual vec3 get_s0_at_scan_point(std::size_t index) const = 0; virtual Probe get_probe() const = 0; + virtual std::string get_probe_name() const = 0; virtual void set_direction(vec3 direction) = 0; virtual void set_wavelength(double wavelength) = 0; @@ -304,6 +306,21 @@ namespace dxtbx { namespace model { return probe_; } + std::string get_probe_name() const { + // Return a name that matches NeXus definitions from + // https://manual.nexusformat.org/classes/base_classes/NXsource.html + switch (probe_) { + case xray: + return std::string("x-ray"); + case electron: + return std::string("electron"); + case neutron: + return std::string("neutron"); + default: + return std::string("unknown"); + } + } + void set_probe(Probe probe) { probe_ = probe; } @@ -412,6 +429,7 @@ namespace dxtbx { namespace model { /** Print beam information */ inline std::ostream &operator<<(std::ostream &os, const Beam &b) { os << "Beam:\n"; + os << " probe: " << b.get_probe_name() << "\n"; os << " wavelength: " << b.get_wavelength() << "\n"; os << " sample to source direction : " << b.get_sample_to_source_direction().const_ref() << "\n"; diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index 370a22b28..6892813e9 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -261,6 +261,7 @@ namespace dxtbx { namespace model { namespace boost_python { .def("get_s0_at_scan_points", &BeamBase::get_s0_at_scan_points) .def("get_s0_at_scan_point", &BeamBase::get_s0_at_scan_point) .def("get_probe", &BeamBase::get_probe) + .def("get_probe_name", &BeamBase::get_probe_name) .def("set_probe", &BeamBase::set_probe) .def("reset_scan_points", &BeamBase::reset_scan_points) .def("rotate_around_origin", From 463391b6964f7a494854e26cef349ab9fa0646d4 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 13:46:23 +0100 Subject: [PATCH 10/26] Include the probe in to/from dict methods and pickle suite. --- src/dxtbx/model/beam.h | 6 ++++-- src/dxtbx/model/beam.py | 4 ++++ src/dxtbx/model/boost_python/beam.cc | 15 +++++++++++---- 3 files changed, 19 insertions(+), 6 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index f3f1ed6a2..b34b31ba0 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -174,6 +174,7 @@ namespace dxtbx { namespace model { * @param polarization_fraction The polarization fraction * @param flux The beam flux * @param transmission The beam transmission + * @param probe The probe value */ Beam(vec3 direction, double wavelength, @@ -182,7 +183,8 @@ namespace dxtbx { namespace model { vec3 polarization_normal, double polarization_fraction, double flux, - double transmission) + double transmission, + Probe probe) : wavelength_(wavelength), divergence_(divergence), sigma_divergence_(sigma_divergence), @@ -190,7 +192,7 @@ namespace dxtbx { namespace model { polarization_fraction_(polarization_fraction), flux_(flux), transmission_(transmission), - probe_(Probe::xray) { + probe_(probe) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index f128e85a1..376b8f0ae 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -101,6 +101,10 @@ def from_dict(d, t=None): joint.update(d) # Create the model from the joint dictionary + if "probe" in joint: + joint["probe"] = Probe.values[joint["probe"]] + else: + joint["probe"] = Probe.xray return Beam.from_dict(joint) @staticmethod diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index 6892813e9..8ad5f6964 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -41,7 +41,8 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get_polarization_normal(), obj.get_polarization_fraction(), obj.get_flux(), - obj.get_transmission()); + obj.get_transmission(), + obj.get_probe()); } static boost::python::tuple getstate(boost::python::object obj) { @@ -108,6 +109,7 @@ namespace dxtbx { namespace model { namespace boost_python { double polarization_fraction, double flux, double transmission, + Probe probe, bool deg) { Beam *beam = NULL; if (deg) { @@ -118,7 +120,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } else { beam = new Beam(sample_to_source, wavelength, @@ -127,7 +130,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } return beam; } @@ -189,6 +193,7 @@ namespace dxtbx { namespace model { namespace boost_python { result["polarization_fraction"] = obj.get_polarization_fraction(); result["flux"] = obj.get_flux(); result["transmission"] = obj.get_transmission(); + result["probe"] = obj.get_probe(); if (obj.get_num_scan_points() > 0) { boost::python::list l; scitbx::af::shared > s0_at_scan_points = obj.get_s0_at_scan_points(); @@ -213,7 +218,8 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get("polarization_normal", vec3(0.0, 1.0, 0.0))), boost::python::extract(obj.get("polarization_fraction", 0.999)), boost::python::extract(obj.get("flux", 0)), - boost::python::extract(obj.get("transmission", 1))); + boost::python::extract(obj.get("transmission", 1)), + boost::python::extract(obj.get("probe", Probe::xray))); if (obj.has_key("s0_at_scan_points")) { boost::python::list s0_at_scan_points = boost::python::extract(obj["s0_at_scan_points"]); @@ -307,6 +313,7 @@ namespace dxtbx { namespace model { namespace boost_python { arg("polarization_fraction"), arg("flux"), arg("transmission"), + arg("probe") = Probe::xray, arg("deg") = true))) .def("__str__", &beam_to_string) .def("to_dict", &to_dict) From 0b08c5db08f4a485d94dd07acb751b08f02049fb Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 14:09:47 +0100 Subject: [PATCH 11/26] Add probe to make_polarized_beam, and fix longstanding bug in that method --- src/dxtbx/model/beam.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 376b8f0ae..2436af348 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -153,6 +153,7 @@ def make_polarized_beam( sigma_divergence=None, flux=None, transmission=None, + probe=Probe.xray, ): assert polarization assert 0.0 <= polarization_fraction <= 1.0 @@ -177,6 +178,7 @@ def make_polarized_beam( float(polarization_fraction), float(flux), float(transmission), + probe, ) elif unit_s0: assert wavelength @@ -189,17 +191,23 @@ def make_polarized_beam( float(polarization_fraction), float(flux), float(transmission), + probe, ) else: assert s0 + sum_sq_s0 = s0[0] ** 2 + s0[1] ** 2 + s0[2] ** 2 + assert sum_sq_s0 > 0 + wavelength = 1.0 / math.sqrt(sum_sq_s0) return Beam( tuple(map(float, s0)), + wavelength, float(divergence), float(sigma_divergence), tuple(map(float, polarization)), float(polarization_fraction), float(flux), float(transmission), + probe, ) @staticmethod From f5e5441fbd1333ea57c4558de90f96f9fd121162 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 14:15:31 +0100 Subject: [PATCH 12/26] Set electron probe for various 3D ED format classes --- src/dxtbx/format/FormatGatanDM4.py | 2 ++ src/dxtbx/format/FormatMRC.py | 2 ++ src/dxtbx/format/FormatSMVTimePix_SU.py | 2 ++ 3 files changed, 6 insertions(+) diff --git a/src/dxtbx/format/FormatGatanDM4.py b/src/dxtbx/format/FormatGatanDM4.py index d5cb31181..acb93f1ef 100644 --- a/src/dxtbx/format/FormatGatanDM4.py +++ b/src/dxtbx/format/FormatGatanDM4.py @@ -21,6 +21,7 @@ ) from dxtbx.format.Format import Format from dxtbx.format.FormatMultiImage import FormatMultiImage +from dxtbx.model.beam import Probe def read_tag(f, byteorder): @@ -358,6 +359,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) def _scan(self): diff --git a/src/dxtbx/format/FormatMRC.py b/src/dxtbx/format/FormatMRC.py index 0b9b06226..b7d9b30bb 100644 --- a/src/dxtbx/format/FormatMRC.py +++ b/src/dxtbx/format/FormatMRC.py @@ -19,6 +19,7 @@ from dxtbx.format.Format import Format from dxtbx.format.FormatMultiImage import FormatMultiImage from dxtbx.model import ScanFactory +from dxtbx.model.beam import Probe logger = logging.getLogger("dials") @@ -226,6 +227,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) diff --git a/src/dxtbx/format/FormatSMVTimePix_SU.py b/src/dxtbx/format/FormatSMVTimePix_SU.py index 223061836..b3492ff56 100644 --- a/src/dxtbx/format/FormatSMVTimePix_SU.py +++ b/src/dxtbx/format/FormatSMVTimePix_SU.py @@ -12,6 +12,7 @@ from scitbx import matrix from dxtbx.format.FormatSMV import FormatSMV +from dxtbx.model.beam import Probe from dxtbx.model.detector import Detector @@ -76,6 +77,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) def _scan(self): From 9f2400ad2cb1fc796e16e4f4ecf7982c6dbc4e3e Mon Sep 17 00:00:00 2001 From: David Waterman Date: Mon, 24 Jul 2023 16:24:38 +0100 Subject: [PATCH 13/26] News --- newsfragments/xxx.feature | 1 + 1 file changed, 1 insertion(+) create mode 100644 newsfragments/xxx.feature diff --git a/newsfragments/xxx.feature b/newsfragments/xxx.feature new file mode 100644 index 000000000..0ff6b0846 --- /dev/null +++ b/newsfragments/xxx.feature @@ -0,0 +1 @@ +The ``Beam`` model now has a ``probe`` value to keep track of the type of radiation. From 823b075aac0c9d5192b06904feefcb3365462a92 Mon Sep 17 00:00:00 2001 From: DiamondLightSource-build-server Date: Mon, 24 Jul 2023 15:32:56 +0000 Subject: [PATCH 14/26] Rename newsfragments/xxx.feature to newsfragments/647.feature --- newsfragments/{xxx.feature => 647.feature} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename newsfragments/{xxx.feature => 647.feature} (100%) diff --git a/newsfragments/xxx.feature b/newsfragments/647.feature similarity index 100% rename from newsfragments/xxx.feature rename to newsfragments/647.feature From 930be0957804d84be170e921e34c9633b92b5e4f Mon Sep 17 00:00:00 2001 From: David Waterman Date: Tue, 25 Jul 2023 09:54:50 +0100 Subject: [PATCH 15/26] Easier override of panel geometry (#644) When setting panels by PHIL, merge multiple definitions first by id. This fixes #299. --- newsfragments/299.bugfix | 1 + src/dxtbx/model/detector.py | 48 +++++++++++++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 2 deletions(-) create mode 100644 newsfragments/299.bugfix diff --git a/newsfragments/299.bugfix b/newsfragments/299.bugfix new file mode 100644 index 000000000..5a8abfc77 --- /dev/null +++ b/newsfragments/299.bugfix @@ -0,0 +1 @@ +Panel geometry definitions in PHIL are merged by panel id before constructing panels diff --git a/src/dxtbx/model/detector.py b/src/dxtbx/model/detector.py index c98828c72..c046ca8d9 100644 --- a/src/dxtbx/model/detector.py +++ b/src/dxtbx/model/detector.py @@ -1,6 +1,7 @@ from __future__ import annotations import os +from copy import deepcopy import pycbf @@ -212,6 +213,42 @@ ) +def merge_panel_scope_extracts_by_id(panel_params): + + id_to_params = {} + for (i, params) in enumerate(panel_params): + if params.id not in id_to_params: + id_to_params[params.id] = [ + i, + ] + else: + id_to_params[params.id].append(i) + + merged_params = [] + for params_set in id_to_params.values(): + params0 = deepcopy(panel_params[params_set[0]]) + for i in params_set[1:]: + params1 = deepcopy(panel_params[i]) + for key in params1.__dict__: + if ( + key.startswith("_") + or key == "id" + or params0.__dict__[key] == params1.__dict__[key] + ): + continue + if ( + params0.__dict__[key] is not None + and params1.__dict__[key] is not None + ): + raise RuntimeError( + f"Multiple definitions for {key} for panel id={params0.id}" + ) + if params0.__dict__[key] is None: + params0.__dict__[key] = params1.__dict__[key] + merged_params.append(params0) + return merged_params + + class DetectorFactory: """A factory class for detector objects, which will encapsulate standard detector designs to make it a little easier to get started with these. In @@ -229,7 +266,11 @@ def generate_from_phil(params, beam=None): # Create a list of panels panel_list = {} - for panel_params in params.detector.panel: + + # merge panel params by id first + merged = merge_panel_scope_extracts_by_id(params.detector.panel) + + for panel_params in merged: panel = Panel() if panel_params.name is not None: panel.set_name(panel_params.name) @@ -335,8 +376,11 @@ def overwrite_from_phil(params, detector, beam=None): Overwrite from phil parameters """ + # merge panel params by id first + merged = merge_panel_scope_extracts_by_id(params.detector.panel) + # Override any panel parameters - for panel_params in params.detector.panel: + for panel_params in merged: panel = detector[panel_params.id] if panel_params.name is not None: panel.set_name(panel_params.name) From 436a1b6a03298338a4027c4ca374b73480adff14 Mon Sep 17 00:00:00 2001 From: Nicholas Devenish Date: Tue, 25 Jul 2023 18:42:06 +0100 Subject: [PATCH 16/26] Flumpy: Don't allow converting 1D vectors to vec2/3 (#648) Previously, a size (3,) array would have erroneously converted to a 0x3 length vec3. This is now explicitly forbidden; you must pass at least a (0, 3) multidimensional array. Fixes #439. --- newsfragments/439.bugfix | 1 + src/dxtbx/boost_python/flumpy.cc | 7 +++++++ tests/test_flumpy.py | 10 +++++++++- 3 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 newsfragments/439.bugfix diff --git a/newsfragments/439.bugfix b/newsfragments/439.bugfix new file mode 100644 index 000000000..a97cc5116 --- /dev/null +++ b/newsfragments/439.bugfix @@ -0,0 +1 @@ +``flumpy``: Fix case where incorrect ``flex.vec2``, ``flex.vec3`` could be generated. diff --git a/src/dxtbx/boost_python/flumpy.cc b/src/dxtbx/boost_python/flumpy.cc index 12e5ebf94..d7ce34b20 100644 --- a/src/dxtbx/boost_python/flumpy.cc +++ b/src/dxtbx/boost_python/flumpy.cc @@ -520,6 +520,13 @@ py::object vec_from_numpy(py::array np_array) { static_assert(VecType::fixed_size == 2 || VecType::fixed_size == 3, "Only vec2/vec3 supported"); + + // Only accept arrays that have a dimension higher than 1 - we want + // numpy.array([1,2,3]) to fail but numpy.array([[1,2,3]]) to work + if (np_array.ndim() == 1) { + throw std::invalid_argument("Array for conversion to vec must be multidimensional"); + } + // Only accept arrays whose last dimension is the size of this object if (np_array.shape(np_array.ndim() - 1) != VecType::fixed_size) { throw std::invalid_argument("Input array last dimension is not size " diff --git a/tests/test_flumpy.py b/tests/test_flumpy.py index cd09e3c09..88ef9f8cc 100644 --- a/tests/test_flumpy.py +++ b/tests/test_flumpy.py @@ -136,7 +136,7 @@ def test_reverse_numeric_2d(flex_numeric): def test_numeric_4d(flex_numeric): - #  Check that we can think fourth-dimnesionally + # Check that we can think fourth-dimensionally grid = flex.grid(1, 9, 8, 5) fo = flex_numeric(grid) assert fo.nd() == 4 @@ -402,3 +402,11 @@ def test_int_long_degeneracy(): npo[0] = 42 assert all(fo[x] == npo[x] for x in range(4)) assert fo[0] == 42 + + +def test_single_entry_vec_from_numpy(): + with pytest.raises(ValueError): + flumpy.vec_from_numpy(np.array([2.0, 3.0, 4.0])) + ao = flumpy.vec_from_numpy(np.array([[2.0, 3.0, 4.0]])) + assert len(ao) == 1 + assert ao[0] == pytest.approx([2.0, 3.0, 4.0]) From 180fac839980dd912abf918be34c99f86d56ab62 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Wed, 26 Jul 2023 14:17:57 +0100 Subject: [PATCH 17/26] FormatROD (#645) Partial support for Rigaku Oxford Diffraction format images. --------- Co-authored-by: Takanori Nakane Co-authored-by: Nicholas Devenish --- newsfragments/645.feature | 1 + src/dxtbx/format/FormatROD.py | 459 ++++++++++++++++++++++++++++++++++ 2 files changed, 460 insertions(+) create mode 100644 newsfragments/645.feature create mode 100644 src/dxtbx/format/FormatROD.py diff --git a/newsfragments/645.feature b/newsfragments/645.feature new file mode 100644 index 000000000..fb38a5121 --- /dev/null +++ b/newsfragments/645.feature @@ -0,0 +1 @@ +Add partial support for the Rigaku Oxford Diffraction file format. diff --git a/src/dxtbx/format/FormatROD.py b/src/dxtbx/format/FormatROD.py new file mode 100644 index 000000000..417f99f1a --- /dev/null +++ b/src/dxtbx/format/FormatROD.py @@ -0,0 +1,459 @@ +"""Support for the Rigaku Oxford Diffraction image format + +NB: Rigaku datasets may use a non-zero-padded image incremental serial number. +At present, this is not compatible with assumptions in dxtbx. In order to +import these datasets, the images should be renumbered first. + +See https://github.com/cctbx/dxtbx/issues/646 for details.""" + +# Takanori Nakane took David Waterman's code to parse headers from +# https://github.com/cctbx/cctbx_project/commit/b95467f3b2a70a37eeb820ea294128a32551700c +# and heavily modified it. The original commit in the cctbx_project repository is orphan now. +from __future__ import annotations + +__author__ = "David Waterman, Takanori Nakane" +__copyright__ = ( + "Copyright 2018 United Kingdom Research and Innovation & 2022 Takanori Nakane" +) +__license__ = "BSD 3-clause" + +import re +import struct + +import numpy as np + +from scitbx.array_family import flex + +from dxtbx.format.Format import Format + + +class FormatROD(Format): + """An image reading class for Rigaku Oxford Diffraction images. + + The Rigaku Oxford Diffraction format is used by CrysAlis(Pro). + """ + + def __init__(self, image_file, **kwargs): + """Initialise the image structure from the given file.""" + + from dxtbx import IncorrectFormatError + + if not self.understand(image_file): + raise IncorrectFormatError(self, image_file) + + Format.__init__(self, image_file, **kwargs) + + @staticmethod + def understand(image_file): + with FormatROD.open_file(image_file, "rb") as f: + try: + hdr = f.read(256).decode("ascii") + except UnicodeDecodeError: + return False + + lines = hdr.splitlines() + if len(lines) < 2: + return False + + vers = lines[0].split() + if len(vers) < 2 or vers[0] != "OD" or vers[1] != "SAPPHIRE": + return False + + compression = lines[1].split("=") + if compression[0] != "COMPRESSION": + return False + + return True + + @staticmethod + def _read_ascii_header(image_file): + """Read the ASCII header comprising the first 256 bytes of the file""" + + hd = {} + with FormatROD.open_file(image_file, "rb") as f: + hdr = f.read(256).decode("ascii") + lines = hdr.splitlines() + + vers = lines[0].split() + if len(vers) < 2 or vers[0] != "OD" or vers[1] != "SAPPHIRE": + raise ValueError("Wrong header format") + hd["version"] = float(vers[-1]) + + compression = lines[1].split("=") + if compression[0] != "COMPRESSION": + raise ValueError("Wrong header format") + hd["compression"] = compression[1] + + # Extract definitions from the 3rd - 5th line + defn = re.compile(r"([A-Z]+=[ 0-9]+)") + for line in lines[2:5]: + sizes = defn.findall(line) + for s in sizes: + n, v = s.split("=") + hd[n] = int(v) + + hd["time"] = lines[5].split("TIME=")[-1].strip("\x1a").rstrip() + + return hd + + @staticmethod + def _read_binary_header( + image_file, + offset=256, + nbytes=5120, + general_nbytes=512, + special_nbytes=768, + km4gonio_nbytes=1024, + statistics_nbytes=512, + history_nbytes=2048, + ): + """Read the most relevant parameters from the binary header""" + + with FormatROD.open_file(image_file, "rb") as f: + # General section + f.seek(offset) + bin_x, bin_y = struct.unpack("> 4) & 15) + ipos += 1 + # ipos_bit = ipos * 8 + + for i in range(2): + nbit = nbits[i] + # 0: 0 + # 1: 1 or 0 + # 2: -1, 0, 1, 2 + # 3: -3, -2, -1, 0, 1, 2, 3, 4 (= b111 -3 = 7 - 3) + zero_at = 0 + if nbit > 1: + zero_at = (1 << (nbit - 1)) - 1 + + v = 0 + for j in range(nbit): + v |= linedata[ipos] << (8 * j) + ipos += 1 + + mask = (1 << nbit) - 1 + for j in range(BLOCKSIZE): + ret[opos] = ((v >> (nbit * j)) & mask) - zero_at + opos += 1 + + for i in range(opos - BLOCKSIZE * 2, opos): + offset = ret[i] + + # surprisingly, this IF is very slow! + if offset >= SHORT_OVERFLOW_SIGNED: + if offset >= LONG_OVERFLOW_SIGNED: + offset = linedata[ipos : (ipos + 4)].view(np.int32)[0] + ipos += 4 + else: + offset = linedata[ipos : (ipos + 2)].view(np.int16)[0] + ipos += 2 + + ret[i] = offset + ret[i - 1] + + for i in range(nrest): + px = linedata[ipos] + ipos += 1 + if px < SHORT_OVERFLOW: + ret[opos] = ret[opos - 1] + px - 127 + elif px == LONG_OVERFLOW: + ret[opos] = ( + ret[opos - 1] + linedata[ipos : (ipos + 4)].view(np.int32)[0] + ) + ipos += 4 + else: + ret[opos] = ( + ret[opos - 1] + linedata[ipos : (ipos + 2)].view(np.int16)[0] + ) + ipos += 2 + opos += 1 + + return ret + + +if __name__ == "__main__": + import sys + + for arg in sys.argv[1:]: + if FormatROD.understand(arg): + reader = FormatROD(arg) + print(reader._txt_header) + print(reader._bin_header) + else: + print("Unsupported format") From e9b88d273389b5fdaa2afc37fa7940905fbb5cb9 Mon Sep 17 00:00:00 2001 From: Nicholas Devenish Date: Wed, 26 Jul 2023 14:41:57 +0100 Subject: [PATCH 18/26] Check Registry against refl/expt files (#650) This was previously done implicitly, but only when dials_regression was present. This makes the test explicit, and works with dials-data (so that it will work on e.g. CI with restricted access to test cases). --- newsfragments/650.misc | 1 + tests/format/test_format.py | 12 ++++++++++++ 2 files changed, 13 insertions(+) create mode 100644 newsfragments/650.misc create mode 100644 tests/format/test_format.py diff --git a/newsfragments/650.misc b/newsfragments/650.misc new file mode 100644 index 000000000..94f0ed982 --- /dev/null +++ b/newsfragments/650.misc @@ -0,0 +1 @@ +Format classes are now tested against invalid binary data with dials-data, for when dials-regression is not present. diff --git a/tests/format/test_format.py b/tests/format/test_format.py new file mode 100644 index 000000000..3b951ba8d --- /dev/null +++ b/tests/format/test_format.py @@ -0,0 +1,12 @@ +from __future__ import annotations + +from dxtbx.format import Registry + + +def test_reading_refl_failure(dials_data): + test_data = dials_data("centroid_test_data", pathlib=True) + + # Without dials_regression, none of the dials-data tests check for this "invalid binary data" case + assert Registry.get_format_class_for_file(test_data / "indexed.refl") is None + # Check .expt while here + assert Registry.get_format_class_for_file(test_data / "indexed.expt") is None From 69169ff5a0e144c445f29677bed0b18fd6d8d0d4 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Wed, 26 Jul 2023 16:58:34 +0100 Subject: [PATCH 19/26] Make Imageset slices consistently index from 0 Without this, slices sometimes needed to be indexed by their original image index, even inside the slice. The previous workaround in 67aab2 (#485) for negative offset was reacting to the erroneous behaviour introduced by including batch offsets in imagesequence slicing. --- newsfragments/633.bugfix | 1 + src/dxtbx/imageset.py | 23 +++-------------------- src/dxtbx/model/experiment_list.py | 3 +-- tests/test_imageset.py | 8 +++++++- 4 files changed, 12 insertions(+), 23 deletions(-) create mode 100644 newsfragments/633.bugfix diff --git a/newsfragments/633.bugfix b/newsfragments/633.bugfix new file mode 100644 index 000000000..7ce5110e5 --- /dev/null +++ b/newsfragments/633.bugfix @@ -0,0 +1 @@ +Slicing of imageset objects is now consistently 0-based, including for the sliced data accessor. Previously, the data accessor had to be accessed with the original index offsets. diff --git a/src/dxtbx/imageset.py b/src/dxtbx/imageset.py index 574779956..85fd2200b 100644 --- a/src/dxtbx/imageset.py +++ b/src/dxtbx/imageset.py @@ -300,29 +300,12 @@ def __getitem__(self, item): if not isinstance(item, slice): return self.get_corrected_data(item) else: - offset = self.get_scan().get_batch_offset() if item.step is not None: raise IndexError("Sequences must be sequential") - # nasty workaround for https://github.com/dials/dials/issues/1153 - # slices with -1 in them are meaningful :-/ so grab the original - # constructor arguments of the slice object. - # item.start and item.stop may have been compromised at this point. - if offset < 0: - start, stop, step = item.__reduce__()[1] - if start is None: - start = 0 - else: - start -= offset - if stop is None: - stop = len(self) - else: - stop -= offset - else: - start = item.start or 0 - stop = item.stop or (len(self) + offset) - start -= offset - stop -= offset + start = item.start or 0 + stop = item.stop or len(self) + if self.data().has_single_file_reader(): reader = self.reader().copy(self.reader().paths(), stop - start) else: diff --git a/src/dxtbx/model/experiment_list.py b/src/dxtbx/model/experiment_list.py index d415e1ccd..ec836614f 100644 --- a/src/dxtbx/model/experiment_list.py +++ b/src/dxtbx/model/experiment_list.py @@ -625,8 +625,7 @@ def from_sequence_and_crystal(imageset, crystal, load_models=True): # if imagesequence is still images, make one experiment for each # all referencing into the same image set if imageset.get_scan().is_still(): - start, end = imageset.get_scan().get_array_range() - for j in range(start, end): + for j in range(len(imageset)): subset = imageset[j : j + 1] experiments.append( Experiment( diff --git a/tests/test_imageset.py b/tests/test_imageset.py index 9285db95d..7c54d65b7 100644 --- a/tests/test_imageset.py +++ b/tests/test_imageset.py @@ -387,6 +387,11 @@ def tst_get_item(self, sequence): with pytest.raises(RuntimeError): _ = sequence2[5] + # Check data access matches expected images from the slice + panel_data1 = sequence[3][0] + panel_data2 = sequence2[0][0] + assert panel_data1.all_eq(panel_data2) + assert len(sequence2) == 4 assert_can_get_detectorbase(sequence2, range(0, 4), 5) self.tst_get_models(sequence2, range(0, 4), 5) @@ -396,11 +401,12 @@ def tst_get_item(self, sequence): with pytest.raises(IndexError): _ = sequence[3:7:2] + # Batch offset should not affect slicing of imagesequence # Simulate a scan starting from image 0 sequence_ = copy.deepcopy(sequence) sequence_.get_scan().set_batch_offset(-1) sequence3 = sequence_[3:7] - assert sequence3.get_array_range() == (4, 8) + assert sequence3.get_array_range() == (3, 7) @staticmethod def tst_paths(sequence, filenames1): From 0a068b56df0cc05cd5ea5e5ebe2f0c203e341e54 Mon Sep 17 00:00:00 2001 From: Nicholas Devenish Date: Thu, 27 Jul 2023 14:52:33 +0100 Subject: [PATCH 20/26] MNT: Update CMake modules path relative to script Without this, dxtbx did not work as a subdirectory in a larger CMakeLists project. --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c286d39e2..96c6b76d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) project(dxtbx LANGUAGES C CXX) # Add the included modules -set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules/") # General cmake environment configuration include(SetDefaultBuildRelWithDebInfo) # Default builds to release with debug info From 64bcb3240a3e300fe7c7fa2f0de7a7aa931d9ef3 Mon Sep 17 00:00:00 2001 From: Christian Becke Date: Mon, 31 Jul 2023 10:58:00 +0200 Subject: [PATCH 21/26] Fixes for EMBL beamlines at PETRA (#626) The E-32-0017 Eiger moved to P13 on 2021-5-22, where the goniometer axis no longer needs overriding. Check for that date in the MiniCBF, and check for the serial matching the date range. Co-authored-by: Christian Becke Co-authored-by: Nicholas Devenish --- newsfragments/626.feature | 1 + src/dxtbx/format/FormatCBFMini.py | 15 +++++++++++++++ src/dxtbx/format/FormatCBFMiniEigerPetraP14.py | 13 ++++++++++++- 3 files changed, 28 insertions(+), 1 deletion(-) create mode 100644 newsfragments/626.feature diff --git a/newsfragments/626.feature b/newsfragments/626.feature new file mode 100644 index 000000000..c28cc4293 --- /dev/null +++ b/newsfragments/626.feature @@ -0,0 +1 @@ +Update Format handling to reflect move of Eiger detector from PETRA P14 to P13. diff --git a/src/dxtbx/format/FormatCBFMini.py b/src/dxtbx/format/FormatCBFMini.py index 4561665fd..39060f84a 100644 --- a/src/dxtbx/format/FormatCBFMini.py +++ b/src/dxtbx/format/FormatCBFMini.py @@ -8,6 +8,7 @@ from __future__ import annotations import binascii +import datetime import os import pathlib import sys @@ -73,6 +74,20 @@ def __init__(self, image_file, **kwargs): self._raw_data = None super().__init__(image_file, **kwargs) + @staticmethod + def _get_timestamp_from_raw_header( + header: str | list[str], + ) -> datetime.datetime | None: + """Given a raw header, or lines from, attempt to extract the timestamp field""" + if isinstance(header, str): + header = header.splitlines() + timestamp = None + for record in header: + if len(record[1:].split()) <= 2 and record.count(":") == 2: + timestamp = datetime.datetime.fromisoformat(record[1:].strip()) + break + return timestamp + def _start(self): """Open the image file, read the image header, copy it into a dictionary for future reference.""" diff --git a/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py b/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py index 530948590..36b5cce91 100644 --- a/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py +++ b/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py @@ -3,6 +3,7 @@ from __future__ import annotations +import datetime import sys from dxtbx.format.FormatCBFMiniEiger import FormatCBFMiniEiger @@ -19,11 +20,21 @@ def understand(image_file): header = FormatCBFMiniEiger.get_cbf_header(image_file) + # Valid from 22nd May 2021 + expected_serial = "E-32-0129" + if timestamp := FormatCBFMiniEiger._get_timestamp_from_raw_header(header): + # We have a timestamp. Let's see what detector we should expect + + # Before 22nd May 2021 + if timestamp < datetime.datetime(2021, 5, 22): + expected_serial = "E-32-0107" + + # Find the line recording detector serial, and check for record in header.split("\n"): if ( "# detector" in record.lower() and "eiger" in record.lower() - and "E-32-0107" in record + and expected_serial in record ): return True From ef5f9abab7c08d46f2bff43ad497e44fe3604048 Mon Sep 17 00:00:00 2001 From: David McDonagh <60879630+toastisme@users.noreply.github.com> Date: Wed, 2 Aug 2023 14:15:44 +0100 Subject: [PATCH 22/26] Add polychromatic beam (#621) Add new Beam class "PolychromaticBeam" for polychromatic/multi-wavelength/wide bandpass experiments. --------- Co-authored-by: DiamondLightSource-build-server Co-authored-by: Nicholas Devenish --- newsfragments/621.feature | 1 + src/dxtbx/dxtbx_model_ext.pyi | 29 ++++ src/dxtbx/model/__init__.py | 3 + src/dxtbx/model/beam.h | 191 +++++++++++++++++++++++++-- src/dxtbx/model/beam.py | 145 +++++++++++++------- src/dxtbx/model/boost_python/beam.cc | 130 ++++++++++++++++++ tests/model/test_beam.py | 94 ++++++++++++- 7 files changed, 535 insertions(+), 58 deletions(-) create mode 100644 newsfragments/621.feature diff --git a/newsfragments/621.feature b/newsfragments/621.feature new file mode 100644 index 000000000..394ff8e74 --- /dev/null +++ b/newsfragments/621.feature @@ -0,0 +1 @@ +Add new Beam class "PolychromaticBeam" for polychromatic/multi-wavelength/wide bandpass experiments. diff --git a/src/dxtbx/dxtbx_model_ext.pyi b/src/dxtbx/dxtbx_model_ext.pyi index 4cbee00fb..4eb5e1530 100644 --- a/src/dxtbx/dxtbx_model_ext.pyi +++ b/src/dxtbx/dxtbx_model_ext.pyi @@ -114,6 +114,35 @@ class Beam(BeamBase): def from_dict(data: Dict) -> Beam: ... def to_dict(self) -> Dict: ... +class PolychromaticBeam(Beam): + @overload + def __init__(self, beam: PolychromaticBeam) -> None: ... + @overload + def __init__(self, direction: Vec3Float) -> None: ... + @overload + def __init__( + self, + direction: Vec3Float, + divergence: float, + sigma_divergence: float, + deg: bool = ..., + ) -> None: ... + @overload + def __init__( + self, + direction: Vec3Float, + divergence: float, + sigma_divergence: float, + polarization_normal: Vec3Float, + polarization_fraction: float, + flux: float, + transmission: float, + deg: bool = ..., + ) -> None: ... + @staticmethod + def from_dict(data: Dict) -> PolychromaticBeam: ... + def to_dict(self) -> Dict: ... + class CrystalBase: @property def num_scan_points(self) -> int: ... diff --git a/src/dxtbx/model/__init__.py b/src/dxtbx/model/__init__.py index 26a2fc4d8..a7e3ada72 100644 --- a/src/dxtbx/model/__init__.py +++ b/src/dxtbx/model/__init__.py @@ -45,6 +45,7 @@ OffsetPxMmStrategy, Panel, ParallaxCorrectedPxMmStrategy, + PolychromaticBeam, PxMmStrategy, Scan, ScanBase, @@ -80,6 +81,7 @@ OffsetPxMmStrategy, Panel, ParallaxCorrectedPxMmStrategy, + PolychromaticBeam, PxMmStrategy, Scan, ScanBase, @@ -97,6 +99,7 @@ __all__ = ( "Beam", "BeamBase", + "PolychromaticBeam", "BeamFactory", "Crystal", "CrystalBase", diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 9903341ec..a8d72f33c 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -188,7 +188,7 @@ namespace dxtbx { namespace model { return direction_; } - double get_wavelength() const { + virtual double get_wavelength() const { return wavelength_; } @@ -207,16 +207,16 @@ namespace dxtbx { namespace model { direction_ = direction.normalize(); } - void set_wavelength(double wavelength) { + virtual void set_wavelength(double wavelength) { wavelength_ = wavelength; } - vec3 get_s0() const { + virtual vec3 get_s0() const { DXTBX_ASSERT(wavelength_ != 0.0); return -direction_ * 1.0 / wavelength_; } - void set_s0(vec3 s0) { + virtual void set_s0(vec3 s0) { DXTBX_ASSERT(s0.length() > 0); direction_ = -s0.normalize(); wavelength_ = 1.0 / s0.length(); @@ -293,7 +293,7 @@ namespace dxtbx { namespace model { s0_at_scan_points_.clear(); } - bool operator==(const BeamBase &rhs) const { + virtual bool operator==(const BeamBase &rhs) const { double eps = 1.0e-6; // scan-varying model checks @@ -327,11 +327,11 @@ namespace dxtbx { namespace model { <= eps; } - bool is_similar_to(const BeamBase &rhs, - double wavelength_tolerance, - double direction_tolerance, - double polarization_normal_tolerance, - double polarization_fraction_tolerance) const { + virtual bool is_similar_to(const BeamBase &rhs, + double wavelength_tolerance, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { // scan varying model checks if (get_num_scan_points() != rhs.get_num_scan_points()) { return false; @@ -375,8 +375,7 @@ namespace dxtbx { namespace model { friend std::ostream &operator<<(std::ostream &os, const Beam &b); - private: - double wavelength_; + protected: vec3 direction_; double divergence_; double sigma_divergence_; @@ -384,6 +383,9 @@ namespace dxtbx { namespace model { double polarization_fraction_; double flux_; double transmission_; + + private: + double wavelength_; scitbx::af::shared > s0_at_scan_points_; }; @@ -403,6 +405,171 @@ namespace dxtbx { namespace model { return os; } + class PolychromaticBeam : public Beam { + public: + PolychromaticBeam() { + set_direction(vec3(0.0, 0.0, 1.0)); + set_divergence(0.0); + set_sigma_divergence(0.0); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + } + + /** + * @param direction The beam direction pointing sample to source + */ + PolychromaticBeam(vec3 direction) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(0.0); + set_sigma_divergence(0.0); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + } + + /** + * @param direction The beam direction pointing sample to source + * @param divergence The beam divergence + * @param sigma_divergence The standard deviation of the beam divergence + */ + PolychromaticBeam(vec3 direction, + double divergence, + double sigma_divergence) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(divergence); + set_sigma_divergence(sigma_divergence); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + } + + /** + * @param direction The beam direction pointing sample to source + * @param divergence The beam divergence + * @param sigma_divergence The standard deviation of the beam divergence + * @param polarization_normal The polarization plane + * @param polarization_fraction The polarization fraction + * @param flux The beam flux + * @param transmission The beam transmission + */ + PolychromaticBeam(vec3 direction, + double divergence, + double sigma_divergence, + vec3 polarization_normal, + double polarization_fraction, + double flux, + double transmission) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(divergence); + set_sigma_divergence(sigma_divergence); + set_polarization_normal(polarization_normal); + set_polarization_fraction(polarization_fraction); + set_flux(flux); + set_transmission(transmission); + } + + double get_wavelength() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed wavelength"); + return -1.; + } + + void set_wavelength(double wavelength) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed wavelength"); + } + + vec3 get_s0() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return vec3(0., 0., 0.); + } + + void set_s0(vec3 s0) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + std::size_t get_num_scan_points() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return 1; + } + + void set_s0_at_scan_points(const scitbx::af::const_ref > &s0) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + scitbx::af::shared > get_s0_at_scan_points() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return scitbx::af::shared >(1, (0., 0., 0.)); + } + + vec3 get_s0_at_scan_point(std::size_t index) const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return vec3(0., 0., 0.); + } + + void reset_scan_points() { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + bool operator==(const BeamBase &rhs) const { + double eps = 1.0e-6; + + return std::abs(angle_safe(direction_, rhs.get_sample_to_source_direction())) + <= eps + && std::abs(divergence_ - rhs.get_divergence()) <= eps + && std::abs(sigma_divergence_ - rhs.get_sigma_divergence()) <= eps + && std::abs( + angle_safe(polarization_normal_, rhs.get_polarization_normal())) + <= eps + && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) + <= eps; + } + + bool is_similar_to(const BeamBase &rhs, + double wavelength_tolerance, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { + return is_similar_to(rhs, + direction_tolerance, + polarization_normal_tolerance, + polarization_fraction_tolerance); + } + + bool is_similar_to(const BeamBase &rhs, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { + return std::abs(angle_safe(direction_, rhs.get_sample_to_source_direction())) + <= direction_tolerance + && std::abs( + angle_safe(polarization_normal_, rhs.get_polarization_normal())) + <= polarization_normal_tolerance + && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) + <= polarization_fraction_tolerance; + } + }; + + /** Print beam information */ + inline std::ostream &operator<<(std::ostream &os, const PolychromaticBeam &b) { + os << "Beam:\n"; + os << " sample to source direction : " + << b.get_sample_to_source_direction().const_ref() << "\n"; + os << " divergence: " << b.get_divergence() << "\n"; + os << " sigma divergence: " << b.get_sigma_divergence() << "\n"; + os << " polarization normal: " << b.get_polarization_normal().const_ref() + << "\n"; + os << " polarization fraction: " << b.get_polarization_fraction() << "\n"; + os << " flux: " << b.get_flux() << "\n"; + os << " transmission: " << b.get_transmission() << "\n"; + return os; + } + }} // namespace dxtbx::model #endif // DXTBX_MODEL_BEAM_H diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 870dd9c4b..7c94cffa7 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -1,15 +1,18 @@ from __future__ import annotations import math +from typing import Tuple import pycbf import libtbx.phil try: - from ..dxtbx_model_ext import Beam + from ..dxtbx_model_ext import Beam, PolychromaticBeam except ModuleNotFoundError: - from dxtbx_model_ext import Beam # type: ignore + from dxtbx_model_ext import Beam, PolychromaticBeam # type: ignore + +Vec3Float = Tuple[float, float, float] beam_phil_scope = libtbx.phil.parse( """ @@ -17,6 +20,10 @@ .expert_level = 1 .short_caption = "Beam overrides" { + type = *monochromatic polychromatic + .type = choice + .help = "Override the beam type" + .short_caption = "beam_type" wavelength = None .type = float @@ -27,6 +34,14 @@ .help = "Override the sample to source direction" .short_caption = "Sample to source direction" + divergence = None + .type = float + .help = "Override the beam divergence" + + sigma_divergence = None + .type = float + .help = "Override the beam sigma divergence" + polarization_normal = None .type = floats(size=3) .help = "Override the polarization normal" @@ -57,62 +72,77 @@ class BeamFactory: will be used, otherwise simplified descriptions can be applied.""" @staticmethod - def from_phil(params, reference=None): + def from_phil( + params: libtbx.phil.scope_extract, + reference: Beam | PolychromaticBeam = None, + ) -> Beam | PolychromaticBeam: """ Convert the phil parameters into a beam model - """ + # Check the input if reference is None: - beam = Beam() + beam = ( + PolychromaticBeam() if params.beam.type == "polychromatic" else Beam() + ) else: beam = reference # Set the parameters - if params.beam.wavelength is not None: - beam.set_wavelength(params.beam.wavelength) - elif reference is None: - raise RuntimeError("No wavelength set") + if params.beam.type == "monochromatic": + if params.beam.wavelength is not None: + beam.set_wavelength(params.beam.wavelength) + elif reference is None: + raise RuntimeError("No wavelength set") if params.beam.direction is not None: beam.set_direction(params.beam.direction) elif reference is None: raise RuntimeError("No beam direction set") + + if params.beam.divergence is not None: + beam.set_divergence(params.beam.divergence) + if params.beam.sigma_divergence is not None: + beam.set_sigma_divergence(params.beam.sigma_divergence) if params.beam.polarization_normal is not None: beam.set_polarization_normal(params.beam.polarization_normal) if params.beam.polarization_fraction is not None: beam.set_polarization_fraction(params.beam.polarization_fraction) + if params.beam.transmission is not None: + beam.set_transmission(params.beam.transmission) + if params.beam.flux is not None: + beam.set_flux(params.beam.flux) return beam @staticmethod - def from_dict(d, t=None): - """Convert the dictionary to a beam model + def from_dict(dict: dict, template: dict = None) -> Beam | PolychromaticBeam: + """Convert the dictionary to a beam model""" - Params: - d The dictionary of parameters - t The template dictionary to use + if template is not None: + if "__id__" in dict and "__id__" in template: + assert ( + dict["__id__"] == template["__id__"] + ), "Beam and template dictionaries are not the same type." - Returns: - The beam model - """ - if d is None and t is None: + if dict is None and template is None: return None - joint = t.copy() if t else {} - joint.update(d) + joint = template.copy() if template else {} + joint.update(dict) # Create the model from the joint dictionary + if joint.get("__id__") == "polychromatic": + return PolychromaticBeam.from_dict(joint) return Beam.from_dict(joint) @staticmethod def make_beam( - sample_to_source=None, - wavelength=None, - s0=None, - unit_s0=None, - divergence=None, - sigma_divergence=None, - ): - + sample_to_source: Vec3Float = None, + wavelength: float = None, + s0: Vec3Float = None, + unit_s0: Vec3Float = None, + divergence: float = None, + sigma_divergence: float = None, + ) -> Beam: if divergence is None or sigma_divergence is None: divergence = 0.0 sigma_divergence = 0.0 @@ -137,19 +167,41 @@ def make_beam( assert s0 return Beam(tuple(map(float, s0))) + @staticmethod + def make_polychromatic_beam( + direction: Vec3Float, + divergence: float = 0.0, + sigma_divergence: float = 0.0, + polarization_normal: Vec3Float = (0.0, 1.0, 0.0), + polarization_fraction: float = 0.5, + flux: float = 0.0, + transmission: float = 1.0, + deg: bool = True, + ) -> PolychromaticBeam: + return PolychromaticBeam( + tuple(map(float, direction)), + float(divergence), + float(sigma_divergence), + tuple(map(float, polarization_normal)), + float(polarization_fraction), + float(flux), + float(transmission), + bool(deg), + ) + @staticmethod def make_polarized_beam( - sample_to_source=None, - wavelength=None, - s0=None, - unit_s0=None, - polarization=None, - polarization_fraction=None, - divergence=None, - sigma_divergence=None, - flux=None, - transmission=None, - ): + sample_to_source: Vec3Float = None, + wavelength: float = None, + s0: Vec3Float = None, + unit_s0: Vec3Float = None, + polarization: Vec3Float = None, + polarization_fraction: float = None, + divergence: float = None, + sigma_divergence: float = None, + flux: float = None, + transmission: float = None, + ) -> Beam: assert polarization assert 0.0 <= polarization_fraction <= 1.0 @@ -199,7 +251,7 @@ def make_polarized_beam( ) @staticmethod - def simple(wavelength): + def simple(wavelength: float) -> Beam: """Construct a beam object on the principle that the beam is aligned with the +z axis, as is quite normal. Also assume the beam has polarization fraction 0.999 and is polarized in the x-z plane, unless @@ -219,7 +271,7 @@ def simple(wavelength): ) @staticmethod - def simple_directional(sample_to_source, wavelength): + def simple_directional(sample_to_source: Vec3Float, wavelength: float) -> Beam: """Construct a beam with direction and wavelength.""" if wavelength > 0.05: @@ -236,8 +288,11 @@ def simple_directional(sample_to_source, wavelength): @staticmethod def complex( - sample_to_source, polarization_fraction, polarization_plane_normal, wavelength - ): + sample_to_source: Vec3Float, + polarization_fraction: float, + polarization_plane_normal: Vec3Float, + wavelength: float, + ) -> Beam: """Full access to the constructor for cases where we do know everything that we need...""" @@ -249,7 +304,7 @@ def complex( ) @staticmethod - def imgCIF(cif_file): + def imgCIF(cif_file: str) -> Beam: """Initialize a detector model from an imgCIF file. N.B. the definition of the polarization plane is not completely helpful in this - it is the angle between the polarization plane and the @@ -263,7 +318,7 @@ def imgCIF(cif_file): return result @staticmethod - def imgCIF_H(cbf_handle): + def imgCIF_H(cbf_handle: pycbf.cbf_handle_struct) -> Beam: """Initialize a detector model from an imgCIF file. N.B. the definition of the polarization plane is not completely helpful in this - it is the angle between the polarization plane and the diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index bf709d526..45bec9f95 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -25,6 +25,8 @@ namespace dxtbx { namespace model { namespace boost_python { using scitbx::deg_as_rad; using scitbx::rad_as_deg; + // Beam + std::string beam_to_string(const Beam &beam) { std::stringstream ss; ss << beam; @@ -180,6 +182,7 @@ namespace dxtbx { namespace model { namespace boost_python { template <> boost::python::dict to_dict(const Beam &obj) { boost::python::dict result; + result["__id__"] = "monochromatic"; result["direction"] = obj.get_sample_to_source_direction(); result["wavelength"] = obj.get_wavelength(); result["divergence"] = rad_as_deg(obj.get_divergence()); @@ -221,6 +224,102 @@ namespace dxtbx { namespace model { namespace boost_python { return b; } + // PolychromaticBeam + + std::string PolychromaticBeam_to_string(const PolychromaticBeam &beam) { + std::stringstream ss; + ss << beam; + return ss.str(); + } + + struct PolychromaticBeamPickleSuite : boost::python::pickle_suite { + static boost::python::tuple getinitargs(const PolychromaticBeam &obj) { + return boost::python::make_tuple(obj.get_sample_to_source_direction(), + obj.get_divergence(), + obj.get_sigma_divergence(), + obj.get_polarization_normal(), + obj.get_polarization_fraction(), + obj.get_flux(), + obj.get_transmission()); + } + }; + + static PolychromaticBeam *make_PolychromaticBeam(vec3 direction) { + return new PolychromaticBeam(direction); + } + + static PolychromaticBeam *make_PolychromaticBeam_w_divergence(vec3 direction, + double divergence, + double sigma_divergence, + bool deg) { + PolychromaticBeam *beam = NULL; + if (deg) { + beam = new PolychromaticBeam( + direction, deg_as_rad(divergence), deg_as_rad(sigma_divergence)); + } else { + beam = new PolychromaticBeam(direction, divergence, sigma_divergence); + } + return beam; + } + + static PolychromaticBeam *make_PolychromaticBeam_w_all( + vec3 direction, + double divergence, + double sigma_divergence, + vec3 polarization_normal, + double polarization_fraction, + double flux, + double transmission, + bool deg) { + PolychromaticBeam *beam = NULL; + if (deg) { + beam = new PolychromaticBeam(direction, + deg_as_rad(divergence), + deg_as_rad(sigma_divergence), + polarization_normal, + polarization_fraction, + flux, + transmission); + } else { + beam = new PolychromaticBeam(direction, + divergence, + sigma_divergence, + polarization_normal, + polarization_fraction, + flux, + transmission); + } + return beam; + } + + template <> + boost::python::dict to_dict(const PolychromaticBeam &obj) { + boost::python::dict result; + result["__id__"] = "polychromatic"; + result["direction"] = obj.get_sample_to_source_direction(); + result["divergence"] = rad_as_deg(obj.get_divergence()); + result["sigma_divergence"] = rad_as_deg(obj.get_sigma_divergence()); + result["polarization_normal"] = obj.get_polarization_normal(); + result["polarization_fraction"] = obj.get_polarization_fraction(); + result["flux"] = obj.get_flux(); + result["transmission"] = obj.get_transmission(); + return result; + } + + template <> + PolychromaticBeam *from_dict(boost::python::dict obj) { + PolychromaticBeam *b = new PolychromaticBeam( + boost::python::extract >(obj["direction"]), + deg_as_rad(boost::python::extract(obj.get("divergence", 0.0))), + deg_as_rad(boost::python::extract(obj.get("sigma_divergence", 0.0))), + boost::python::extract >( + obj.get("polarization_normal", vec3(0.0, 1.0, 0.0))), + boost::python::extract(obj.get("polarization_fraction", 0.999)), + boost::python::extract(obj.get("flux", 0)), + boost::python::extract(obj.get("transmission", 1))); + return b; + } + void export_beam() { using namespace beam_detail; @@ -305,6 +404,37 @@ namespace dxtbx { namespace model { namespace boost_python { .staticmethod("from_dict") .def_pickle(BeamPickleSuite()); + class_, bases >( + "PolychromaticBeam") + .def("__init__", + make_constructor( + &make_PolychromaticBeam, default_call_policies(), (arg("direction")))) + .def("__init__", + make_constructor(&make_PolychromaticBeam_w_divergence, + default_call_policies(), + (arg("direction"), + arg("divergence"), + arg("sigma_divergence"), + arg("deg") = true))) + .def("__init__", + make_constructor(&make_PolychromaticBeam_w_all, + default_call_policies(), + (arg("direction"), + arg("divergence"), + arg("sigma_divergence"), + arg("polarization_normal"), + arg("polarization_fraction"), + arg("flux"), + arg("transmission"), + arg("deg") = true))) + .def("__str__", &PolychromaticBeam_to_string) + .def("to_dict", &to_dict) + .def("from_dict", + &from_dict, + return_value_policy()) + .staticmethod("from_dict") + .def_pickle(PolychromaticBeamPickleSuite()); + scitbx::af::boost_python::flex_wrapper::plain("flex_Beam"); } diff --git a/tests/model/test_beam.py b/tests/model/test_beam.py index 7f1ef5b9d..2f83a790c 100644 --- a/tests/model/test_beam.py +++ b/tests/model/test_beam.py @@ -5,7 +5,7 @@ from libtbx.phil import parse from scitbx import matrix -from dxtbx.model import Beam +from dxtbx.model import Beam, PolychromaticBeam from dxtbx.model.beam import BeamFactory, beam_phil_scope @@ -176,3 +176,95 @@ def test_beam_object_comparison(): def test_beam_self_serialization(): beam = Beam() assert beam == BeamFactory.from_dict(beam.to_dict()) + + +def test_polychromatic_beam_from_phil(): + params = beam_phil_scope.fetch( + parse( + """ + beam { + type = polychromatic + direction = (0., 0., 1.) + divergence = 0.2 + sigma_divergence = 0.3 + polarization_normal = (0., -1., 0.) + polarization_fraction = .65 + transmission = .5 + flux = .75 + } + """ + ) + ).extract() + + beam = BeamFactory.from_phil(params) + assert isinstance(beam, PolychromaticBeam) + + assert beam.get_sample_to_source_direction() == pytest.approx((0.0, 0.0, 1.0)) + assert beam.get_divergence() == pytest.approx(0.2) + assert beam.get_sigma_divergence() == pytest.approx(0.3) + assert beam.get_polarization_normal() == pytest.approx((0.0, -1.0, 0.0)) + assert beam.get_polarization_fraction() == pytest.approx(0.65) + assert beam.get_transmission() == pytest.approx(0.5) + assert beam.get_flux() == pytest.approx(0.75) + + +def test_polychromatic_beam_from_dict(): + beam = PolychromaticBeam() + assert beam == BeamFactory.from_dict(beam.to_dict()) + + +def test_make_polychromatic_beam(): + + direction = (0.0, 0.0, 1.0) + divergence = 0.2 + sigma_divergence = 0.3 + polarization_normal = (0.0, -1.0, 0.0) + polarization_fraction = 0.65 + transmission = 0.5 + flux = 0.75 + + beam = BeamFactory.make_polychromatic_beam( + direction=direction, + divergence=divergence, + sigma_divergence=sigma_divergence, + polarization_normal=polarization_normal, + polarization_fraction=polarization_fraction, + transmission=transmission, + flux=flux, + ) + + assert beam.get_sample_to_source_direction() == pytest.approx((0.0, 0.0, 1.0)) + assert beam.get_divergence() == pytest.approx(0.2) + assert beam.get_sigma_divergence() == pytest.approx(0.3) + assert beam.get_polarization_normal() == pytest.approx((0.0, -1.0, 0.0)) + assert beam.get_polarization_fraction() == pytest.approx(0.65) + assert beam.get_transmission() == pytest.approx(0.5) + assert beam.get_flux() == pytest.approx(0.75) + + +def test_polychromatic_beam_wavelength_guards(): + beam = PolychromaticBeam() + with pytest.raises(RuntimeError): + _ = beam.get_wavelength() + with pytest.raises(RuntimeError): + _ = beam.get_s0() + with pytest.raises(RuntimeError): + _ = beam.get_num_scan_points() + with pytest.raises(RuntimeError): + _ = beam.get_s0_at_scan_points() + with pytest.raises(RuntimeError): + _ = beam.get_s0_at_scan_point(0) + with pytest.raises(RuntimeError): + beam.reset_scan_points() + with pytest.raises(RuntimeError): + beam.set_wavelength(1.0) + with pytest.raises(RuntimeError): + beam.set_s0((0.0, 0.0, 0.1)) + + +def test_polychromatic_beam_str(): + beam = PolychromaticBeam() + assert ( + beam.__str__() + == "Beam:\n sample to source direction : {0,0,1}\n divergence: 0\n sigma divergence: 0\n polarization normal: {0,1,0}\n polarization fraction: 0.5\n flux: 0\n transmission: 1\n" + ) From d2cf031c5ceaebb218d15acf7ec363fe86a53327 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Wed, 2 Aug 2023 17:11:14 +0100 Subject: [PATCH 23/26] Add probe to PolychromaticBeam --- src/dxtbx/model/beam.h | 17 +++++++++++++---- src/dxtbx/model/boost_python/beam.cc | 15 +++++++++++---- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index b6083f592..799c7f45c 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -423,11 +423,11 @@ namespace dxtbx { namespace model { double polarization_fraction_; double flux_; double transmission_; + Probe probe_; private: double wavelength_; scitbx::af::shared > s0_at_scan_points_; - Probe probe_; }; /** Print beam information */ @@ -457,6 +457,7 @@ namespace dxtbx { namespace model { set_polarization_fraction(0.5); set_flux(0); set_transmission(1.0); + set_probe(Probe::xray); } /** @@ -471,6 +472,7 @@ namespace dxtbx { namespace model { set_polarization_fraction(0.5); set_flux(0); set_transmission(1.0); + set_probe(Probe::xray); } /** @@ -489,6 +491,7 @@ namespace dxtbx { namespace model { set_polarization_fraction(0.5); set_flux(0); set_transmission(1.0); + set_probe(Probe::xray); } /** @@ -499,6 +502,7 @@ namespace dxtbx { namespace model { * @param polarization_fraction The polarization fraction * @param flux The beam flux * @param transmission The beam transmission + * @param probe The probe value */ PolychromaticBeam(vec3 direction, double divergence, @@ -506,7 +510,8 @@ namespace dxtbx { namespace model { vec3 polarization_normal, double polarization_fraction, double flux, - double transmission) { + double transmission, + Probe probe) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); set_divergence(divergence); @@ -515,6 +520,7 @@ namespace dxtbx { namespace model { set_polarization_fraction(polarization_fraction); set_flux(flux); set_transmission(transmission); + set_probe(probe); } double get_wavelength() const { @@ -569,7 +575,8 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= eps && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= eps; + <= eps + && (probe_ == rhs.get_probe()); } bool is_similar_to(const BeamBase &rhs, @@ -593,13 +600,15 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= polarization_normal_tolerance && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= polarization_fraction_tolerance; + <= polarization_fraction_tolerance + && (probe_ == rhs.get_probe()); } }; /** Print beam information */ inline std::ostream &operator<<(std::ostream &os, const PolychromaticBeam &b) { os << "Beam:\n"; + os << " probe: " << b.get_probe_name() << "\n"; os << " sample to source direction : " << b.get_sample_to_source_direction().const_ref() << "\n"; os << " divergence: " << b.get_divergence() << "\n"; diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index 173fc96d4..98ca40957 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -247,7 +247,8 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get_polarization_normal(), obj.get_polarization_fraction(), obj.get_flux(), - obj.get_transmission()); + obj.get_transmission(), + obj.get_probe()); } }; @@ -277,6 +278,7 @@ namespace dxtbx { namespace model { namespace boost_python { double polarization_fraction, double flux, double transmission, + Probe probe, bool deg) { PolychromaticBeam *beam = NULL; if (deg) { @@ -286,7 +288,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } else { beam = new PolychromaticBeam(direction, divergence, @@ -294,7 +297,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } return beam; } @@ -310,6 +314,7 @@ namespace dxtbx { namespace model { namespace boost_python { result["polarization_fraction"] = obj.get_polarization_fraction(); result["flux"] = obj.get_flux(); result["transmission"] = obj.get_transmission(); + result["probe"] = obj.get_probe(); return result; } @@ -323,7 +328,8 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get("polarization_normal", vec3(0.0, 1.0, 0.0))), boost::python::extract(obj.get("polarization_fraction", 0.999)), boost::python::extract(obj.get("flux", 0)), - boost::python::extract(obj.get("transmission", 1))); + boost::python::extract(obj.get("transmission", 1)), + boost::python::extract(obj.get("probe", Probe::xray))); return b; } @@ -442,6 +448,7 @@ namespace dxtbx { namespace model { namespace boost_python { arg("polarization_fraction"), arg("flux"), arg("transmission"), + arg("probe") = Probe::xray, arg("deg") = true))) .def("__str__", &PolychromaticBeam_to_string) .def("to_dict", &to_dict) From 7f0e28905c6f364f1e7a39649e117c54a7d0bbd7 Mon Sep 17 00:00:00 2001 From: David Waterman Date: Wed, 2 Aug 2023 17:43:17 +0100 Subject: [PATCH 24/26] Fixes for tests --- src/dxtbx/model/beam.py | 6 ++++-- tests/model/test_beam.py | 7 +++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 10741b433..5ae08b3d5 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -8,9 +8,9 @@ import libtbx.phil try: - from ..dxtbx_model_ext import Beam, Probe, PolychromaticBeam + from ..dxtbx_model_ext import Beam, PolychromaticBeam, Probe except ModuleNotFoundError: - from dxtbx_model_ext import Beam, Probe, PolychromaticBeam # type: ignore + from dxtbx_model_ext import Beam, PolychromaticBeam, Probe # type: ignore Vec3Float = Tuple[float, float, float] @@ -181,6 +181,7 @@ def make_polychromatic_beam( polarization_fraction: float = 0.5, flux: float = 0.0, transmission: float = 1.0, + probe=Probe.xray, deg: bool = True, ) -> PolychromaticBeam: return PolychromaticBeam( @@ -191,6 +192,7 @@ def make_polychromatic_beam( float(polarization_fraction), float(flux), float(transmission), + probe, bool(deg), ) diff --git a/tests/model/test_beam.py b/tests/model/test_beam.py index 2f83a790c..48a6eb5a1 100644 --- a/tests/model/test_beam.py +++ b/tests/model/test_beam.py @@ -6,7 +6,7 @@ from scitbx import matrix from dxtbx.model import Beam, PolychromaticBeam -from dxtbx.model.beam import BeamFactory, beam_phil_scope +from dxtbx.model.beam import BeamFactory, Probe, beam_phil_scope def test_setting_direction_and_wavelength(): @@ -222,6 +222,7 @@ def test_make_polychromatic_beam(): polarization_fraction = 0.65 transmission = 0.5 flux = 0.75 + probe = Probe.neutron beam = BeamFactory.make_polychromatic_beam( direction=direction, @@ -231,6 +232,7 @@ def test_make_polychromatic_beam(): polarization_fraction=polarization_fraction, transmission=transmission, flux=flux, + probe=probe, ) assert beam.get_sample_to_source_direction() == pytest.approx((0.0, 0.0, 1.0)) @@ -240,6 +242,7 @@ def test_make_polychromatic_beam(): assert beam.get_polarization_fraction() == pytest.approx(0.65) assert beam.get_transmission() == pytest.approx(0.5) assert beam.get_flux() == pytest.approx(0.75) + assert beam.get_probe() == Probe.neutron def test_polychromatic_beam_wavelength_guards(): @@ -266,5 +269,5 @@ def test_polychromatic_beam_str(): beam = PolychromaticBeam() assert ( beam.__str__() - == "Beam:\n sample to source direction : {0,0,1}\n divergence: 0\n sigma divergence: 0\n polarization normal: {0,1,0}\n polarization fraction: 0.5\n flux: 0\n transmission: 1\n" + == "Beam:\n probe: x-ray\n sample to source direction : {0,0,1}\n divergence: 0\n sigma divergence: 0\n polarization normal: {0,1,0}\n polarization fraction: 0.5\n flux: 0\n transmission: 1\n" ) From 1f3c257296796dd3c4f1432c58c914a0d5078d6e Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Thu, 3 Aug 2023 13:01:20 +0100 Subject: [PATCH 25/26] Added get_probe_from_name to Beam. Added Probe to beam_phil_scope. Changed Probe in .expt files to be human readable. --- src/dxtbx/model/beam.h | 16 ++++++++++++++++ src/dxtbx/model/beam.py | 10 ++++++++-- src/dxtbx/model/boost_python/beam.cc | 12 ++++++++---- tests/model/test_beam.py | 12 ++++++++++++ 4 files changed, 44 insertions(+), 6 deletions(-) diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 799c7f45c..fe2fda864 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -323,6 +323,22 @@ namespace dxtbx { namespace model { } } + static Probe get_probe_from_name(const std::string probe) { + // Return a Probe matched to NeXus definitions from + // https://manual.nexusformat.org/classes/base_classes/NXsource.html + + if (probe == "x-ray") { + return Probe::xray; + } else if (probe == "electron") { + return Probe::electron; + } else if (probe == "neutron") { + return Probe::neutron; + } + + DXTBX_ERROR("Unknown probe " + probe); + return Probe::xray; + } + void set_probe(Probe probe) { probe_ = probe; } diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 5ae08b3d5..4cd19be32 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -25,6 +25,11 @@ .help = "Override the beam type" .short_caption = "beam_type" + probe = *x-ray electron neutron + .type = choice + .help = "Override the beam probe" + .short_caption = "beam_probe" + wavelength = None .type = float .help = "Override the beam wavelength" @@ -111,6 +116,7 @@ def from_phil( beam.set_transmission(params.beam.transmission) if params.beam.flux is not None: beam.set_flux(params.beam.flux) + beam.set_probe(Beam.get_probe_from_name(params.beam.probe)) return beam @@ -131,9 +137,9 @@ def from_dict(dict: dict, template: dict = None) -> Beam | PolychromaticBeam: # Create the model from the joint dictionary if "probe" in joint: - joint["probe"] = Probe.values[joint["probe"]] + joint["probe"] = joint["probe"] else: - joint["probe"] = Probe.xray + joint["probe"] = Beam.get_probe_name(Probe.xray) if joint.get("__id__") == "polychromatic": return PolychromaticBeam.from_dict(joint) diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index 98ca40957..e5024d9d5 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -196,7 +196,7 @@ namespace dxtbx { namespace model { namespace boost_python { result["polarization_fraction"] = obj.get_polarization_fraction(); result["flux"] = obj.get_flux(); result["transmission"] = obj.get_transmission(); - result["probe"] = obj.get_probe(); + result["probe"] = obj.get_probe_name(); if (obj.get_num_scan_points() > 0) { boost::python::list l; scitbx::af::shared > s0_at_scan_points = obj.get_s0_at_scan_points(); @@ -222,7 +222,8 @@ namespace dxtbx { namespace model { namespace boost_python { boost::python::extract(obj.get("polarization_fraction", 0.999)), boost::python::extract(obj.get("flux", 0)), boost::python::extract(obj.get("transmission", 1)), - boost::python::extract(obj.get("probe", Probe::xray))); + Beam::get_probe_from_name( + boost::python::extract(obj.get("probe", "x-ray")))); if (obj.has_key("s0_at_scan_points")) { boost::python::list s0_at_scan_points = boost::python::extract(obj["s0_at_scan_points"]); @@ -314,7 +315,7 @@ namespace dxtbx { namespace model { namespace boost_python { result["polarization_fraction"] = obj.get_polarization_fraction(); result["flux"] = obj.get_flux(); result["transmission"] = obj.get_transmission(); - result["probe"] = obj.get_probe(); + result["probe"] = obj.get_probe_name(); return result; } @@ -329,7 +330,8 @@ namespace dxtbx { namespace model { namespace boost_python { boost::python::extract(obj.get("polarization_fraction", 0.999)), boost::python::extract(obj.get("flux", 0)), boost::python::extract(obj.get("transmission", 1)), - boost::python::extract(obj.get("probe", Probe::xray))); + Beam::get_probe_from_name( + boost::python::extract(obj.get("probe", "x-ray")))); return b; } @@ -424,6 +426,8 @@ namespace dxtbx { namespace model { namespace boost_python { .def("to_dict", &to_dict) .def("from_dict", &from_dict, return_value_policy()) .staticmethod("from_dict") + .def("get_probe_from_name", &Beam::get_probe_from_name) + .staticmethod("get_probe_from_name") .def_pickle(BeamPickleSuite()); class_, bases >( diff --git a/tests/model/test_beam.py b/tests/model/test_beam.py index 48a6eb5a1..40659f4d4 100644 --- a/tests/model/test_beam.py +++ b/tests/model/test_beam.py @@ -102,6 +102,18 @@ def test_from_phil(): assert b3.get_polarization_fraction() == 0.5 assert b3.get_polarization_normal() == (1.0, 0.0, 0.0) + params3 = beam_phil_scope.fetch( + parse( + """ + beam { + probe = electron + } + """ + ) + ).extract() + b4 = BeamFactory.from_phil(params3, reference) + assert b4.get_probe() == Probe.electron + def test_scan_varying(): direction = matrix.col((0.013142, 0.002200, 1.450476)) From 4312daf256975efa3f7e359b7d11e00f791afdf1 Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Thu, 3 Aug 2023 14:48:33 +0100 Subject: [PATCH 26/26] Added more typehints for beam.py. Added empty PolychromaticBeam constructor at the Python level. Fixed typo in Beam.from_dict. --- src/dxtbx/dxtbx_model_ext.pyi | 4 ++++ src/dxtbx/model/beam.py | 12 +++++------- src/dxtbx/model/boost_python/beam.cc | 1 + 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/dxtbx/dxtbx_model_ext.pyi b/src/dxtbx/dxtbx_model_ext.pyi index 4eb5e1530..8a4a85934 100644 --- a/src/dxtbx/dxtbx_model_ext.pyi +++ b/src/dxtbx/dxtbx_model_ext.pyi @@ -22,6 +22,8 @@ from scitbx.array_family import shared as flex_shared # Attempt to use the stub typing for flex-inheritance from scitbx.array_family.flex import FlexPlain +from dxtbx_model_ext import Probe # type: ignore + # TypeVar for the set of Experiment models that can be joint-accepted # - profile, imageset and scalingmodel are handled as 'object' TExperimentModel = TypeVar( @@ -113,6 +115,8 @@ class Beam(BeamBase): @staticmethod def from_dict(data: Dict) -> Beam: ... def to_dict(self) -> Dict: ... + @staticmethod + def get_probe_from_name(name: str) -> Probe: ... class PolychromaticBeam(Beam): @overload diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 4cd19be32..a7360e9e5 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -79,7 +79,7 @@ class BeamFactory: @staticmethod def from_phil( params: libtbx.phil.scope_extract, - reference: Beam | PolychromaticBeam = None, + reference: Beam | PolychromaticBeam | None = None, ) -> Beam | PolychromaticBeam: """ Convert the phil parameters into a beam model @@ -136,10 +136,8 @@ def from_dict(dict: dict, template: dict = None) -> Beam | PolychromaticBeam: joint.update(dict) # Create the model from the joint dictionary - if "probe" in joint: - joint["probe"] = joint["probe"] - else: - joint["probe"] = Beam.get_probe_name(Probe.xray) + if "probe" not in joint: + joint["probe"] = "x-ray" if joint.get("__id__") == "polychromatic": return PolychromaticBeam.from_dict(joint) @@ -187,7 +185,7 @@ def make_polychromatic_beam( polarization_fraction: float = 0.5, flux: float = 0.0, transmission: float = 1.0, - probe=Probe.xray, + probe: Probe = Probe.xray, deg: bool = True, ) -> PolychromaticBeam: return PolychromaticBeam( @@ -214,7 +212,7 @@ def make_polarized_beam( sigma_divergence: float = None, flux: float = None, transmission: float = None, - probe=Probe.xray, + probe: Probe = Probe.xray, ) -> Beam: assert polarization assert 0.0 <= polarization_fraction <= 1.0 diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index e5024d9d5..86cdab927 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -432,6 +432,7 @@ namespace dxtbx { namespace model { namespace boost_python { class_, bases >( "PolychromaticBeam") + .def(init()) .def("__init__", make_constructor( &make_PolychromaticBeam, default_call_policies(), (arg("direction"))))