Allow Euler convention setting for calibration

Some users strongly prefer to use a different Euler convention for
calibration, such as extrinsic XYZ. It is easy to convert angles from
one Euler convention to another. However, it is difficult to convert
boundary conditions (min/max ranges for each angle) from one convention
to another.

So we add the ability to choose the Euler convention to use in calibration.
This allows the min/max ranges to be specified in lmfit using that convention.

Signed-off-by: Patrick Avery <[email protected]>

hexrd/fitting/calibration/instrument.py

@@ -6,6 +6,8 @@
 from .lmfit_param_handling import (
     add_engineering_constraints,
     create_instr_params,
+    DEFAULT_EULER_CONVENTION,
+    update_instrument_from_params,
     validate_params_list,
 )
 
@@ -19,7 +21,8 @@ def _normalized_ssqr(resd):
 
 class InstrumentCalibrator:
     def __init__(self, *args, engineering_constraints=None,
-                 set_refinements_from_instrument_flags=True):
+                 set_refinements_from_instrument_flags=True,
+                 euler_convention=DEFAULT_EULER_CONVENTION):
         """
         Model for instrument calibration class as a function of
 
@@ -42,6 +45,7 @@ def __init__(self, *args, engineering_constraints=None,
             assert calib.instr is self.instr, \
                 "all calibrators must refer to the same instrument"
         self._engineering_constraints = engineering_constraints
+        self.euler_convention = euler_convention
 
         self.params = self.make_lmfit_params()
         if set_refinements_from_instrument_flags:
@@ -52,7 +56,10 @@ def __init__(self, *args, engineering_constraints=None,
                                       nan_policy='omit')
 
     def make_lmfit_params(self):
-        params = create_instr_params(self.instr)
+        params = create_instr_params(
+            self.instr,
+            euler_convention=self.euler_convention,
+        )
 
         for calibrator in self.calibrators:
             # We pass the params to the calibrator so it can ensure it
@@ -71,7 +78,12 @@ def make_lmfit_params(self):
 
     def update_all_from_params(self, params):
         # Update instrument and material from the lmfit parameters
-        self.instr.update_from_lmfit_parameter_list(params)
+        update_instrument_from_params(
+            self.instr,
+            params,
+            self.euler_convention,
+        )
+
         for calibrator in self.calibrators:
             calibrator.update_from_lmfit_params(params)
 

hexrd/fitting/calibration/laue.py

@@ -14,6 +14,7 @@
 from .calibrator import Calibrator
 from .lmfit_param_handling import (
     create_grain_params,
+    DEFAULT_EULER_CONVENTION,
     rename_to_avoid_collision,
 )
 
@@ -23,18 +24,16 @@ class LaueCalibrator(Calibrator):
 
     def __init__(self, instr, material, grain_params, default_refinements=None,
                  min_energy=5, max_energy=25, tth_distortion=None,
-                 calibration_picks=None):
+                 calibration_picks=None,
+                 euler_convention=DEFAULT_EULER_CONVENTION):
         self.instr = instr
         self.material = material
         self.grain_params = grain_params
         self.default_refinements = default_refinements
         self.energy_cutoffs = [min_energy, max_energy]
         self.tth_distortion = tth_distortion
 
-        self.euler_convention = {
-            'axes_order': 'zxz',
-            'extrinsic': False,
-        }
+        self.euler_convention = euler_convention
 
         self.data_dict = None
         if calibration_picks is not None:

hexrd/fitting/calibration/lmfit_param_handling.py

@@ -1,11 +1,24 @@
+import lmfit
 import numpy as np
 
-from hexrd.instrument import calc_angles_from_beam_vec
-from hexrd.rotations import RotMatEuler
+from hexrd.instrument import (
+    calc_angles_from_beam_vec,
+    calc_beam_vec,
+)
+from hexrd.rotations import (
+    expMapOfQuat,
+    make_rmat_euler,
+    quatOfRotMat,
+    RotMatEuler,
+)
 from hexrd.material.unitcell import _lpname
 
 
-def create_instr_params(instr):
+# First is the axes_order, second is extrinsic
+DEFAULT_EULER_CONVENTION = ('zxz', False)
+
+
+def create_instr_params(instr, euler_convention=DEFAULT_EULER_CONVENTION):
     # add with tuples: (NAME VALUE VARY MIN  MAX  EXPR  BRUTE_STEP)
     parms_list = []
     if instr.has_multi_beam:
@@ -40,31 +53,18 @@ def create_instr_params(instr):
     parms_list.append(('instr_tvec_x', instr.tvec[0], False, -np.inf, np.inf))
     parms_list.append(('instr_tvec_y', instr.tvec[1], False, -np.inf, np.inf))
     parms_list.append(('instr_tvec_z', instr.tvec[2], False, -np.inf, np.inf))
-    euler_convention = {'axes_order': 'zxz',
-                        'extrinsic': False}
     for det_name, panel in instr.detectors.items():
         det = det_name.replace('-', '_')
-        rmat = panel.rmat
-        rme = RotMatEuler(np.zeros(3,),
-                          **euler_convention)
-        rme.rmat = rmat
-        euler = np.degrees(rme.angles)
-
-        parms_list.append((f'{det}_euler_z',
-                           euler[0],
-                           False,
-                           euler[0]-2,
-                           euler[0]+2))
-        parms_list.append((f'{det}_euler_xp',
-                           euler[1],
-                           False,
-                           euler[1]-2,
-                           euler[1]+2))
-        parms_list.append((f'{det}_euler_zpp',
-                           euler[2],
-                           False,
-                           euler[2]-2,
-                           euler[2]+2))
+
+        angles = detector_angles_euler(panel, euler_convention)
+        angle_names = param_names_euler_convention(det, euler_convention)
+
+        for name, angle in zip(angle_names, angles):
+            parms_list.append((name,
+                               angle,
+                               False,
+                               angle - 2,
+                               angle + 2))
 
         parms_list.append((f'{det}_tvec_x',
                            panel.tvec[0],
@@ -93,6 +93,63 @@ def create_instr_params(instr):
     return parms_list
 
 
+def update_instrument_from_params(instr, params, euler_convention):
+    """
+    this function updates the instrument from the
+    lmfit parameter list. we don't have to keep track
+    of the position numbers as the variables are named
+    variables. this will become the standard in the
+    future since bound constraints can be very easily
+    implemented.
+    """
+    if not isinstance(params, lmfit.Parameters):
+        msg = ('Only lmfit.Parameters is acceptable input. '
+               f'Received: {params}')
+        raise NotImplementedError(msg)
+
+    instr.beam_energy = params['beam_energy'].value
+
+    azim = params['beam_azimuth'].value
+    pola = params['beam_polar'].value
+    instr.beam_vector = calc_beam_vec(azim, pola)
+
+    chi = np.radians(params['instr_chi'].value)
+    instr.chi = chi
+
+    instr_tvec = [params['instr_tvec_x'].value,
+                  params['instr_tvec_y'].value,
+                  params['instr_tvec_z'].value]
+    instr.tvec = np.r_[instr_tvec]
+
+    for det_name, detector in instr.detectors.items():
+        det = det_name.replace('-', '_')
+        set_detector_angles_euler(detector, det, params, euler_convention)
+
+        tvec = np.r_[params[f'{det}_tvec_x'].value,
+                     params[f'{det}_tvec_y'].value,
+                     params[f'{det}_tvec_z'].value]
+        detector.tvec = tvec
+        if detector.detector_type.lower() == 'cylindrical':
+            rad = params[f'{det}_radius'].value
+            detector.radius = rad
+
+        distortion_str = f'{det}_distortion_param'
+        if any(distortion_str in p for p in params):
+            if detector.distortion is None:
+                raise RuntimeError(f"distortion discrepancy for '{det}'!")
+            else:
+                names = np.sort([p for p in params if distortion_str in p])
+                distortion = np.r_[[params[n].value for n in names]]
+                try:
+                    detector.distortion.params = distortion
+                except AssertionError:
+                    raise RuntimeError(
+                        f"distortion for '{det}' "
+                        f"expects {len(detector.distortion.params)} "
+                        f"params but got {len(distortion)}"
+                    )
+
+
 def create_tth_parameters(meas_angles):
     parms_list = []
     for ii, tth in enumerate(meas_angles):
@@ -262,3 +319,66 @@ def validate_params_list(params_list):
     if duplicate_names:
         msg = f'Duplicate names found in params list: {duplicate_names}'
         raise LmfitValidationException(msg)
+
+
+EULER_PARAM_NAMES_MAPPING = {
+    None: ('expmap_x', 'expmap_y', 'expmap_z'),
+    ('xyz', True): ('euler_x', 'euler_y', 'euler_z'),
+    ('zxz', False): ('euler_z', 'euler_xp', 'euler_zpp'),
+}
+
+
+def normalize_euler_convention(euler_convention):
+    if isinstance(euler_convention, dict):
+        return (
+            euler_convention['axes_order'],
+            euler_convention['extrinsic'],
+        )
+
+    return euler_convention
+
+
+def param_names_euler_convention(base, euler_convention):
+    normalized = normalize_euler_convention(euler_convention)
+    return [f'{base}_{x}' for x in EULER_PARAM_NAMES_MAPPING[normalized]]
+
+
+def detector_angles_euler(panel, euler_convention):
+    if euler_convention is None:
+        # Return exponential map parameters
+        return panel.tilt
+
+    normalized = normalize_euler_convention(euler_convention)
+    rmat = panel.rmat
+    rme = RotMatEuler(
+        np.zeros(3,),
+        axes_order=normalized[0],
+        extrinsic=normalized[1],
+    )
+
+    rme.rmat = rmat
+    return np.degrees(rme.angles)
+
+
+def set_detector_angles_euler(panel, base_name, params, euler_convention):
+    normalized = normalize_euler_convention(euler_convention)
+    names = param_names_euler_convention(base_name, euler_convention)
+
+    angles = []
+    for name in names:
+        angles.append(params[name].value)
+
+    angles = np.asarray(angles)
+
+    if euler_convention is None:
+        # No conversion needed
+        panel.tilt = angles
+        return
+
+    rmat = make_rmat_euler(
+        np.radians(angles),
+        axes_order=normalized[0],
+        extrinsic=normalized[1],
+    )
+
+    panel.tilt = expMapOfQuat(quatOfRotMat(rmat))

hexrd/fitting/calibration/structureless.py

@@ -5,6 +5,8 @@
     add_engineering_constraints,
     create_instr_params,
     create_tth_parameters,
+    DEFAULT_EULER_CONVENTION,
+    update_instrument_from_params,
 )
 
 
@@ -31,18 +33,20 @@ def __init__(self,
                  instr,
                  data,
                  tth_distortion=None,
-                 engineering_constraints=None):
+                 engineering_constraints=None,
+                 euler_convention=DEFAULT_EULER_CONVENTION):
 
         self._instr = instr
         self._data = data
         self._tth_distortion = tth_distortion
         self._engineering_constraints = engineering_constraints
+        self.euler_convention = euler_convention
         self.make_lmfit_params()
         self.set_minimizer()
 
     def make_lmfit_params(self):
         params = []
-        params += create_instr_params(self.instr)
+        params += create_instr_params(self.instr, self.euler_convention)
         params += create_tth_parameters(self.meas_angles)
 
         params_dict = lmfit.Parameters()
@@ -53,7 +57,12 @@ def make_lmfit_params(self):
         return params_dict
 
     def calc_residual(self, params):
-        self.instr.update_from_lmfit_parameter_list(params)
+        update_instrument_from_params(
+            self.instr,
+            params,
+            self.euler_convention,
+        )
+
         residual = np.empty([0,])
         for ii, (rng, corr_rng) in enumerate(zip(self.meas_angles,
                                                  self.tth_correction)):