diff --git a/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.cc b/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.cc new file mode 100644 index 0000000000..dedb46f74b --- /dev/null +++ b/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.cc @@ -0,0 +1,114 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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(0.0, 1.0, 0.0); + axes[2] = scitbx::vec3(0.0, 0.0, 1.0); + compose(); +} + +scitbx::af::shared SimpleUParameterisation::get_params() { + return params; +} +scitbx::mat3 SimpleUParameterisation::get_state() { + return U_; +} +void SimpleUParameterisation::set_params(scitbx::af::shared p) { + assert(p.size() == 3); + params = p; + compose(); +} +scitbx::af::shared> 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 B) { + orientation_ = cctbx::crystal_orientation(B, true); +} + +scitbx::af::small SymmetrizeReduceEnlarge::forward_independent_parameters() { + Bconverter.forward(orientation_); + return constraints_.independent_params(Bconverter.G); +} + +cctbx::crystal_orientation SymmetrizeReduceEnlarge::backward_orientation( + scitbx::af::small independent) { + scitbx::sym_mat3 ustar = constraints_.all_params(independent); + Bconverter.validate_and_setG(ustar); + orientation_ = Bconverter.back_as_orientation(); + return orientation_; +} + +scitbx::af::shared> SymmetrizeReduceEnlarge::forward_gradients() { + return dB_dp(Bconverter, constraints_); +} + +void SimpleCellParameterisation::compose() { + scitbx::af::small 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 X = SRE.forward_independent_parameters(); + params = scitbx::af::shared(X.size()); + for (int i = 0; i < X.size(); ++i) { + params[i] = 1E5 * X[i]; + } + compose(); +} + +scitbx::af::shared SimpleCellParameterisation::get_params() { + return params; +} +scitbx::mat3 SimpleCellParameterisation::get_state() { + return B_; +} +void SimpleCellParameterisation::set_params(scitbx::af::shared p) { + params = p; + compose(); +} +scitbx::af::shared> SimpleCellParameterisation::get_dS_dp() { + return dS_dp; +} \ No newline at end of file diff --git a/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.h b/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.h new file mode 100644 index 0000000000..8bd2b6bc75 --- /dev/null +++ b/src/dials/algorithms/profile_model/ellipsoid/UBparameterisation.h @@ -0,0 +1,89 @@ +#ifndef DIALS_ALGORITHMS_UBPARAMETERISATION_H +#define DIALS_ALGORITHMS_UBPARAMETERISATION_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class SimpleUParameterisation { +public: + SimpleUParameterisation(const dxtbx::model::Crystal &crystal); + scitbx::af::shared get_params(); + void set_params(scitbx::af::shared p); + scitbx::mat3 get_state(); + scitbx::af::shared> get_dS_dp(); + +private: + scitbx::af::shared params{3, 0.0}; + scitbx::af::shared> axes{3, scitbx::vec3(1.0, 0.0, 0.0)}; + void compose(); + scitbx::mat3 istate{}; + scitbx::mat3 U_{}; + scitbx::af::shared> dS_dp{ + 3, + scitbx::mat3(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 B); + scitbx::af::small forward_independent_parameters(); + cctbx::crystal_orientation backward_orientation( + scitbx::af::small independent); + scitbx::af::shared> forward_gradients(); + +private: + cctbx::sgtbx::space_group space_group_; + cctbx::sgtbx::tensor_rank_2::constraints constraints_; + cctbx::crystal_orientation orientation_{}; + rstbx::symmetry::AG Bconverter{}; +}; + +class SimpleCellParameterisation { +public: + SimpleCellParameterisation(const dxtbx::model::Crystal &crystal); + scitbx::af::shared get_params(); + void set_params(scitbx::af::shared); + scitbx::mat3 get_state(); + scitbx::af::shared> get_dS_dp(); + +private: + scitbx::af::shared params; + void compose(); + scitbx::mat3 B_{}; + scitbx::af::shared> 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", no_init) + .def(init()) + .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", no_init) + .def(init()) + .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 \ No newline at end of file diff --git a/src/dials/algorithms/profile_model/ellipsoid/test_MLtarget_setup.py b/src/dials/algorithms/profile_model/ellipsoid/test_MLtarget_setup.py new file mode 100644 index 0000000000..6d9ec09f33 --- /dev/null +++ b/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") diff --git a/src/dials/algorithms/profile_model/ellipsoid/test_equivalence_of_conditonal.py b/src/dials/algorithms/profile_model/ellipsoid/test_equivalence_of_conditonal.py new file mode 100644 index 0000000000..f9e977b502 --- /dev/null +++ b/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() diff --git a/src/dials/algorithms/profile_model/ellipsoid/test_fn.py b/src/dials/algorithms/profile_model/ellipsoid/test_fn.py new file mode 100644 index 0000000000..915345d910 --- /dev/null +++ b/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)