diff --git a/hexrd/fitting/calibration/instrument.py b/hexrd/fitting/calibration/instrument.py index ca75bcbdf..7ddc7ad72 100644 --- a/hexrd/fitting/calibration/instrument.py +++ b/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) diff --git a/hexrd/fitting/calibration/laue.py b/hexrd/fitting/calibration/laue.py index c3d1781c8..131d68bfe 100644 --- a/hexrd/fitting/calibration/laue.py +++ b/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,7 +24,8 @@ 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 @@ -31,10 +33,7 @@ def __init__(self, instr, material, grain_params, default_refinements=None, 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: diff --git a/hexrd/fitting/calibration/lmfit_param_handling.py b/hexrd/fitting/calibration/lmfit_param_handling.py index f91d30cbe..e89b98106 100644 --- a/hexrd/fitting/calibration/lmfit_param_handling.py +++ b/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)) diff --git a/hexrd/fitting/calibration/structureless.py b/hexrd/fitting/calibration/structureless.py index cddb13c56..af50b6352 100644 --- a/hexrd/fitting/calibration/structureless.py +++ b/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)): diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py index cf24577f7..621a33726 100644 --- a/hexrd/instrument/hedm_instrument.py +++ b/hexrd/instrument/hedm_instrument.py @@ -42,8 +42,6 @@ import h5py -import lmfit - import numpy as np from io import IOBase @@ -70,13 +68,7 @@ from hexrd import xrdutil from hexrd.material.crystallography import PlaneData from hexrd import constants as ct -from hexrd.rotations import ( - angleAxisOfRotMat, - RotMatEuler, - make_rmat_euler, - expMapOfQuat, - quatOfRotMat - ) +from hexrd.rotations import angleAxisOfRotMat, RotMatEuler from hexrd import distortion as distortion_pkg from hexrd.utils.compatibility import h5py_read_string from hexrd.utils.concurrent import distribute_tasks @@ -1118,70 +1110,6 @@ def update_from_parameter_list(self, p): ii += dpnp return - def update_from_lmfit_parameter_list(self, params): - """ - 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) - - self.beam_energy = params['beam_energy'].value - - azim = params['beam_azimuth'].value - pola = params['beam_polar'].value - self.beam_vector = calc_beam_vec(azim, pola) - - chi = np.radians(params['instr_chi'].value) - self.chi = chi - - instr_tvec = [params['instr_tvec_x'].value, - params['instr_tvec_y'].value, - params['instr_tvec_z'].value] - self.tvec = np.r_[instr_tvec] - - for det_name, detector in self.detectors.items(): - det = det_name.replace('-', '_') - euler = np.r_[params[f'{det}_euler_z'].value, - params[f'{det}_euler_xp'].value, - params[f'{det}_euler_zpp'].value] - - rmat = make_rmat_euler(np.radians(euler), - 'zxz', - extrinsic=False) - tilt = expMapOfQuat(quatOfRotMat(rmat)) - detector.tilt = tilt - - 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 extract_polar_maps(self, plane_data, imgser_dict, active_hkls=None, threshold=None, tth_tol=None, eta_tol=0.25):