Extra cpp files for UB

src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.cc

@@ -0,0 +1,114 @@
+#include <iostream>
+#include <scitbx/vec2.h>
+#include <scitbx/vec3.h>
+#include <dxtbx/model/panel.h>
+#include <cmath>
+#include <dials/algorithms/profile_model/ellipsoid/UBparameterisation.h>
+#include <dials/algorithms/profile_model/ellipsoid/g_gradients.h>
+#include <dials/algorithms/refinement/parameterisation/parameterisation_helpers.h>
+#include <dials/array_family/reflection_table.h>
+#include <dxtbx/model/experiment.h>
+#include <map>
+#include <string>
+#include <cctbx/sgtbx/space_group.h>
+#include <cctbx/sgtbx/tensor_rank_2.h>
+#include <cctbx/crystal_orientation.h>
+#include <rstbx/symmetry/constraints/a_g_conversion.h>
+
+void SimpleUParameterisation::compose() {
+  dials::refinement::CrystalOrientationCompose coc(
+    istate, params[0], axes[0], params[1], axes[1], params[2], axes[2]);
+  U_ = coc.U();
+  dS_dp[0] = coc.dU_dphi1();
+  dS_dp[1] = coc.dU_dphi2();
+  dS_dp[2] = coc.dU_dphi3();
+}
+
+SimpleUParameterisation::SimpleUParameterisation(const dxtbx::model::Crystal &crystal) {
+  istate = crystal.get_U();
+  axes[1] = scitbx::vec3<double>(0.0, 1.0, 0.0);
+  axes[2] = scitbx::vec3<double>(0.0, 0.0, 1.0);
+  compose();
+}
+
+scitbx::af::shared<double> SimpleUParameterisation::get_params() {
+  return params;
+}
+scitbx::mat3<double> SimpleUParameterisation::get_state() {
+  return U_;
+}
+void SimpleUParameterisation::set_params(scitbx::af::shared<double> p) {
+  assert(p.size() == 3);
+  params = p;
+  compose();
+}
+scitbx::af::shared<scitbx::mat3<double>> SimpleUParameterisation::get_dS_dp() {
+  return dS_dp;
+}
+
+SymmetrizeReduceEnlarge::SymmetrizeReduceEnlarge(cctbx::sgtbx::space_group space_group)
+    : space_group_(space_group), constraints_(space_group, true) {}
+
+void SymmetrizeReduceEnlarge::set_orientation(scitbx::mat3<double> B) {
+  orientation_ = cctbx::crystal_orientation(B, true);
+}
+
+scitbx::af::small<double, 6> SymmetrizeReduceEnlarge::forward_independent_parameters() {
+  Bconverter.forward(orientation_);
+  return constraints_.independent_params(Bconverter.G);
+}
+
+cctbx::crystal_orientation SymmetrizeReduceEnlarge::backward_orientation(
+  scitbx::af::small<double, 6> independent) {
+  scitbx::sym_mat3<double> ustar = constraints_.all_params(independent);
+  Bconverter.validate_and_setG(ustar);
+  orientation_ = Bconverter.back_as_orientation();
+  return orientation_;
+}
+
+scitbx::af::shared<scitbx::mat3<double>> SymmetrizeReduceEnlarge::forward_gradients() {
+  return dB_dp(Bconverter, constraints_);
+}
+
+void SimpleCellParameterisation::compose() {
+  scitbx::af::small<double, 6> vals(params.size());
+  for (int i = 0; i < params.size(); ++i) {
+    vals[i] = 1E-5 * params[i];
+  }
+  SRE.set_orientation(B_);
+  SRE.forward_independent_parameters();
+  B_ = SRE.backward_orientation(vals).reciprocal_matrix();
+  dS_dp = SRE.forward_gradients();
+  for (int i = 0; i < dS_dp.size(); ++i) {
+    for (int j = 0; j < 9; ++j) {
+      dS_dp[i][j] *= 1E-5;
+    }
+  }
+}
+
+SimpleCellParameterisation::SimpleCellParameterisation(
+  const dxtbx::model::Crystal &crystal)
+    : B_(crystal.get_B()), SRE(crystal.get_space_group()) {
+  // first get params
+  SRE.set_orientation(B_);
+  scitbx::af::small<double, 6> X = SRE.forward_independent_parameters();
+  params = scitbx::af::shared<double>(X.size());
+  for (int i = 0; i < X.size(); ++i) {
+    params[i] = 1E5 * X[i];
+  }
+  compose();
+}
+
+scitbx::af::shared<double> SimpleCellParameterisation::get_params() {
+  return params;
+}
+scitbx::mat3<double> SimpleCellParameterisation::get_state() {
+  return B_;
+}
+void SimpleCellParameterisation::set_params(scitbx::af::shared<double> p) {
+  params = p;
+  compose();
+}
+scitbx::af::shared<scitbx::mat3<double>> SimpleCellParameterisation::get_dS_dp() {
+  return dS_dp;
+}

src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.h

@@ -0,0 +1,89 @@
+#ifndef DIALS_ALGORITHMS_UBPARAMETERISATION_H
+#define DIALS_ALGORITHMS_UBPARAMETERISATION_H
+
+#include <boost/python.hpp>
+#include <boost/python/def.hpp>
+#include <scitbx/vec2.h>
+#include <scitbx/vec3.h>
+#include <scitbx/mat2.h>
+#include <scitbx/mat3.h>
+#include <cctbx/sgtbx/space_group.h>
+#include <cctbx/sgtbx/tensor_rank_2.h>
+#include <cctbx/crystal_orientation.h>
+#include <rstbx/symmetry/constraints/a_g_conversion.h>
+#include <dxtbx/model/experiment.h>
+
+class SimpleUParameterisation {
+public:
+  SimpleUParameterisation(const dxtbx::model::Crystal &crystal);
+  scitbx::af::shared<double> get_params();
+  void set_params(scitbx::af::shared<double> p);
+  scitbx::mat3<double> get_state();
+  scitbx::af::shared<scitbx::mat3<double>> get_dS_dp();
+
+private:
+  scitbx::af::shared<double> params{3, 0.0};
+  scitbx::af::shared<scitbx::vec3<double>> axes{3, scitbx::vec3<double>(1.0, 0.0, 0.0)};
+  void compose();
+  scitbx::mat3<double> istate{};
+  scitbx::mat3<double> U_{};
+  scitbx::af::shared<scitbx::mat3<double>> dS_dp{
+    3,
+    scitbx::mat3<double>(1.0, 0, 0, 0, 1.0, 0, 0, 0, 1.0)};
+};
+
+class SymmetrizeReduceEnlarge {
+public:
+  SymmetrizeReduceEnlarge(cctbx::sgtbx::space_group space_group);
+  void set_orientation(scitbx::mat3<double> B);
+  scitbx::af::small<double, 6> forward_independent_parameters();
+  cctbx::crystal_orientation backward_orientation(
+    scitbx::af::small<double, 6> independent);
+  scitbx::af::shared<scitbx::mat3<double>> forward_gradients();
+
+private:
+  cctbx::sgtbx::space_group space_group_;
+  cctbx::sgtbx::tensor_rank_2::constraints<double> constraints_;
+  cctbx::crystal_orientation orientation_{};
+  rstbx::symmetry::AG Bconverter{};
+};
+
+class SimpleCellParameterisation {
+public:
+  SimpleCellParameterisation(const dxtbx::model::Crystal &crystal);
+  scitbx::af::shared<double> get_params();
+  void set_params(scitbx::af::shared<double>);
+  scitbx::mat3<double> get_state();
+  scitbx::af::shared<scitbx::mat3<double>> get_dS_dp();
+
+private:
+  scitbx::af::shared<double> params;
+  void compose();
+  scitbx::mat3<double> B_{};
+  scitbx::af::shared<scitbx::mat3<double>> dS_dp{};
+  SymmetrizeReduceEnlarge SRE;
+};
+
+using namespace boost::python;
+
+namespace dials { namespace algorithms { namespace boost_python {
+
+  BOOST_PYTHON_MODULE(dials_algorithms_profile_model_ellipsoid_UBparameterisation_ext) {
+    class_<SimpleUParameterisation>("SimpleUParameterisation", no_init)
+      .def(init<dxtbx::model::Crystal>())
+      .def("get_params", &SimpleUParameterisation::get_params)
+      .def("set_params", &SimpleUParameterisation::set_params)
+      .def("get_state", &SimpleUParameterisation::get_state)
+      .def("get_dS_dp", &SimpleUParameterisation::get_dS_dp);
+
+    class_<SimpleCellParameterisation>("SimpleCellParameterisation", no_init)
+      .def(init<dxtbx::model::Crystal>())
+      .def("get_params", &SimpleCellParameterisation::get_params)
+      .def("set_params", &SimpleCellParameterisation::set_params)
+      .def("get_state", &SimpleCellParameterisation::get_state)
+      .def("get_dS_dp", &SimpleCellParameterisation::get_dS_dp);
+  }
+
+}}}  // namespace dials::algorithms::boost_python
+
+#endif  // DIALS_ALGORITHMS_UBPARAMETERISATION_H

src/dials/algorithms/profile_model/ellipsoid/test_MLtarget_setup.py

@@ -0,0 +1,39 @@
+from __future__ import annotations
+
+from dxtbx.serialize import load
+
+from dials.array_family import flex
+
+# from dials_algorithms_profile_model_ellipsoid_refiner_ext import RefinerData
+from dials_algorithms_profile_model_ellipsoid_parameterisation_ext import (
+    ModelState,
+    Simple6MosaicityParameterisation,
+    WavelengthSpreadParameterisation,
+)
+
+expt = load.experiment_list("../indexed.expt")[0]
+refls = flex.reflection_table.from_file("../indexed.refl").split_by_experiment_id()[0]
+
+
+# rd = RefinerData(expt, refls)
+
+wave = WavelengthSpreadParameterisation(0.01)
+print(wave.sigma())
+print(wave.get_param())
+M = Simple6MosaicityParameterisation(
+    flex.double([0.001, 0.001, 0.002, 0.0001, 0.0001, 0.0001])
+)
+print(hex(id(M)))
+print(M.sigma())
+print(list(M.get_params()))
+model = ModelState(expt.crystal, M, wave, False, False, False, False)
+print(list(model.active_parameters()))
+print(model.mosaicity_covariance_matrix())
+print(model.n_active_parameters())
+new_p = flex.double(
+    [0.011, 0.012, 0.013, 10.5, 0.0015, 0.0016, 0.0017, 0.0002, 0.0003, 0.0004, 0.009]
+)
+print(new_p.size())
+model.set_active_parameters(new_p)
+print(list(model.active_parameters()))
+print("done")

src/dials/algorithms/profile_model/ellipsoid/test_equivalence_of_conditonal.py

@@ -0,0 +1,111 @@
+from __future__ import annotations
+
+import numpy as np
+
+from dxtbx import flumpy
+from scitbx import matrix
+
+from dials.algorithms.profile_model.ellipsoid.refiner import ConditionalDistribution
+from dials.array_family import flex
+
+
+def test_equivalence():
+
+    # copy of some test data.
+    norm_s0 = 0.726686
+    mu = (1.4654e-05, 4.0118e-05, 7.2654e-01)
+    dmu = flex.vec3_double(
+        [
+            (-6.04538430e-06, -4.79266282e-05, 2.58915729e-04),
+            (4.82067272e-05, -5.92986582e-06, 3.20202693e-05),
+            (0.0, 0.0, 0.0),
+            (0.0, 0.0, 0.0),
+        ]
+    )
+    S = matrix.sqr(
+        (
+            1.57194461e-07,
+            -3.03137672e-10,
+            1.64014291e-09,
+            -3.03137672e-10,
+            1.60196971e-07,
+            2.83399363e-08,
+            1.64014291e-09,
+            2.83399363e-08,
+            1.23415987e-08,
+        )
+    )
+    Snp = np.array(
+        [
+            1.57194461e-07,
+            -3.03137672e-10,
+            1.64014291e-09,
+            -3.03137672e-10,
+            1.60196971e-07,
+            2.83399363e-08,
+            1.64014291e-09,
+            2.83399363e-08,
+            1.23415987e-08,
+        ]
+    ).reshape(3, 3)
+
+    dS = flex.mat3_double(
+        [
+            (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
+            (
+                7.92955130e-04,
+                -8.07998127e-07,
+                4.12441448e-06,
+                -8.07998127e-07,
+                8.25616327e-11,
+                -2.00825882e-10,
+                4.12441448e-06,
+                -2.00825882e-10,
+                -1.66681314e-10,
+            ),
+            (
+                -8.15627488e-09,
+                7.21646318e-05,
+                -3.89854794e-04,
+                7.21646318e-05,
+                -2.86566012e-07,
+                1.52621558e-06,
+                -3.89854794e-04,
+                1.52621558e-06,
+                -8.12676428e-06,
+            ),
+        ]
+    )
+
+    c = ConditionalDistribution(
+        norm_s0,
+        np.array(mu).reshape(3, 1),
+        flumpy.to_numpy(dmu).reshape(4, 3).transpose(),
+        Snp,
+        flumpy.to_numpy(dS)
+        .reshape(4, 3, 3)
+        .transpose(1, 2, 0),  # expects 3x3x4, need to do it this way
+    )
+
+    derivs_list = np.array([])
+    for d in c.first_derivatives_of_mean():
+        derivs_list = np.append(derivs_list, d.flatten())
+    sigmas_list = np.array([])
+    for d in c.first_derivatives_of_sigma():
+        sigmas_list = np.append(sigmas_list, d.flatten())
+
+    from dials_algorithms_profile_model_ellipsoid_refiner_ext import test_conditional
+
+    print(c.mean())
+    print(c.sigma())
+    print("dmu")
+    print(derivs_list)
+    print("dsigma")
+    print(sigmas_list)
+
+    test_conditional(norm_s0, mu, dmu, S, dS)
+
+
+if __name__ == "__main__":
+    test_equivalence()

src/dials/algorithms/profile_model/ellipsoid/test_fn.py

@@ -0,0 +1,13 @@
+from __future__ import annotations
+
+from scitbx import matrix
+from scitbx.array_family import flex
+
+from dials.algorithms.profile_model.ellipsoid import rse
+
+x = matrix.sqr([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
+y = (1.0, 1.0)
+z = (2.0, 2.0)
+
+res = rse(x, y, z)
+print(res)