diff --git a/Framework/PythonInterface/plugins/algorithms/AlignComponents.py b/Framework/PythonInterface/plugins/algorithms/AlignComponents.py index 04e4a321266e75307255dd72965849a5e546c8fa..78105372ae2dfabd8cbf337a167840368a5735a1 100644 --- a/Framework/PythonInterface/plugins/algorithms/AlignComponents.py +++ b/Framework/PythonInterface/plugins/algorithms/AlignComponents.py @@ -515,641 +515,6 @@ class AlignComponents(PythonAlgorithm): indexes_and_titles = [(index, title) for index, title in enumerate(table_tofs.getColumnNames()) if '@' in title] column_indexes, titles = list(zip(*indexes_and_titles)) peak_tofs = np.array([table_tofs.column(i) for i in column_indexes]) # shape = (peak_count, detector_count) - # sort by increasing p# Mantid Repository : https://github.com/mantidproject/mantid - # # - # # Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI, - # # NScD Oak Ridge National Laboratory, European Spallation Source, - # # Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS - # # SPDX - License - Identifier: GPL - 3.0 + - # #pylint: disable=no-init, no-name-in-module - # import enum - # import math - # import numpy as np - # from typing import List - # - # from mantid.api import ( - # AlgorithmFactory, FileAction, FileProperty, InstrumentValidator, ITableWorkspaceProperty, - # MatrixWorkspaceProperty, Progress, PropertyMode, PythonAlgorithm, WorkspaceProperty) - # from mantid.dataobjects import MaskWorkspace, TableWorkspace - # from mantid.kernel import ( - # Direction, EnabledWhenProperty, FloatArrayProperty, FloatBoundedValidator, logger, PropertyCriterion, Quat, - # StringArrayProperty, StringListValidator, V3D) - # import mantid.simpleapi as api - # - # - # class AlignComponents(PythonAlgorithm): - # """ - # Class to align components - # """ - # - # _optionsList = ["Xposition", "Yposition", "Zposition", "AlphaRotation", "BetaRotation", "GammaRotation"] - # adjustment_items = ['ComponentName', 'Xposition', 'Yposition', 'Zposition', - # 'XdirectionCosine', 'YdirectionCosine', 'ZdirectionCosine', 'RotationAngle'] - # _optionsDict = {} - # _initialPos = None - # _move = False - # _rotate = False - # _eulerConvention = None - # - # def category(self): - # """ - # Mantid required - # """ - # return "Diffraction" - # - # def seeAlso(self): - # return ["GetDetOffsetsMultiPeaks", "CalibrateRectangularDetectors"] - # - # def name(self): - # """ - # Mantid required - # """ - # return "AlignComponents" - # - # def summary(self): - # """ - # Mantid required - # """ - # return "Align a component by minimising difference to an offset workspace" - # - # #pylint: disable=too-many-locals - # def PyInit(self): - # - # # - # # Reference and input data - # self.declareProperty(ITableWorkspaceProperty( - # "PeakCentersTofTable", "", optional=PropertyMode.Mandatory, direction=Direction.Input), - # doc="Table of found peak centers, in TOF units") - # - # self.declareProperty(FloatArrayProperty( - # "PeakPositions", values=[], direction=Direction.Input), - # doc="Comma separated list of reference peak center d-spacings, sorted by increasing value.") - # - # properties = ["PeakCentersTofTable", "PeakPositions"] - # [self.setPropertyGroup(name, 'Reference and Input Data') for name in properties] - # - # # - # # Selection of the instrument and mask - # self.declareProperty( - # MatrixWorkspaceProperty("MaskWorkspace", "", optional=PropertyMode.Optional, direction=Direction.Input), - # doc="Mask workspace") - # - # self.declareProperty(FileProperty( - # name="InstrumentFilename", defaultValue="", action=FileAction.OptionalLoad, extensions=[".xml"]), - # doc="Instrument filename") - # - # self.declareProperty( - # WorkspaceProperty("InputWorkspace", "", validator=InstrumentValidator(), optional=PropertyMode.Optional, - # direction=Direction.Input), - # doc="Workspace containing the instrument to be calibrated") - # - # self.declareProperty( - # WorkspaceProperty("OutputWorkspace", "", optional=PropertyMode.Mandatory, direction=Direction.Output), - # doc='Workspace containing the calibrated instrument') - # - # properties = ['MaskWorkspace', "InstrumentFilename", "InputWorkspace", "OutputWorkspace"] - # [self.setPropertyGroup(name, 'Instrument Options') for name in properties] - # - # # - # # Components - # self.declareProperty( - # "AdjustmentsTable", "", direction=Direction.Input, - # doc="Name of output table containing optimized locations and orientations for each component") - # - # self.declareProperty( - # name="FitSourcePosition", defaultValue=False, - # doc="Fit the source position, changes L1 (source to sample) distance. " - # "Uses entire instrument. Occurs before Components are Aligned.") - # - # self.declareProperty( - # name="FitSamplePosition", defaultValue=False, - # doc="Fit the sample position, changes L1 (source to sample) and L2 (sample to detector) distance." - # "Uses entire instrument. Occurs before Components are Aligned.") - # - # self.declareProperty( - # StringArrayProperty("ComponentList", direction=Direction.Input), - # doc="Comma separated list on instrument components to refine.") - # - # properties = ["AdjustmentsTable", "FitSourcePosition", "FitSamplePosition", "ComponentList"] - # [self.setPropertyGroup(name, 'Declaration of Components') for name in properties] - # - # # - # # Translation Properties - # # X position - # self.declareProperty(name="Xposition", defaultValue=False, - # doc="Refine Xposition of source and/or sample and/or components") - # condition = EnabledWhenProperty("Xposition", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinXposition", defaultValue=-0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Minimum relative X bound (m)") - # self.setPropertySettings("MinXposition", condition) - # self.declareProperty(name="MaxXposition", defaultValue=0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Maximum relative X bound (m)") - # self.setPropertySettings("MaxXposition", condition) - # - # - # self.declareProperty(name="Yposition", defaultValue=False, - # doc="Refine Yposition of source and/or sample and/or components") - # condition = EnabledWhenProperty("Yposition", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinYposition", defaultValue=-0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Minimum relative Y bound (m)") - # self.setPropertySettings("MinYposition", condition) - # self.declareProperty(name="MaxYposition", defaultValue=0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Maximum relative Y bound (m)") - # self.setPropertySettings("MaxYposition", condition) - # - # # Z position - # self.declareProperty(name="Zposition", defaultValue=False, - # doc="Refine Zposition of source and/or sample and/or components") - # condition = EnabledWhenProperty("Zposition", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinZposition", defaultValue=-0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Minimum relative Z bound (m)") - # self.setPropertySettings("MinZposition", condition) - # self.declareProperty(name="MaxZposition", defaultValue=0.1, - # validator=FloatBoundedValidator(-10.0, 10.0), - # doc="Maximum relative Z bound (m)") - # self.setPropertySettings("MaxZposition", condition) - # - # properties = ["Xposition", "MinXposition", "MaxXposition", "Yposition", "MinYposition", "MaxYposition", - # "Zposition", "MinZposition", "MaxZposition"] - # [self.setPropertyGroup(name, "Translation") for name in properties] - # - # - # # - # # Rotation Properties - # eulerConventions = ["ZXZ", "XYX", "YZY", "ZYZ", "XZX", "YXY", "XYZ", "YZX", "ZXY", "XZY", "ZYX", "YXZ"] - # self.declareProperty( - # name="EulerConvention", defaultValue="YZX", validator=StringListValidator(eulerConventions), - # doc="Euler angles convention used when calculating and displaying angles," - # "eg XYZ corresponding to alpha beta gamma.") - # - # # alpha rotation - # self.declareProperty(name="AlphaRotation", defaultValue=False, - # doc="Refine rotation around first axis, alpha, for the components") - # condition = EnabledWhenProperty("AlphaRotation", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinAlphaRotation", defaultValue=-10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Minimum relative alpha rotation (deg)") - # self.setPropertySettings("MinAlphaRotation", condition) - # self.declareProperty(name="MaxAlphaRotation", defaultValue=10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Maximum relative alpha rotation (deg)") - # self.setPropertySettings("MaxAlphaRotation", condition) - # - # # beta rotation - # self.declareProperty(name="BetaRotation", defaultValue=False, - # doc="Refine rotation around seconds axis, beta, for the components") - # condition = EnabledWhenProperty("BetaRotation", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinBetaRotation", defaultValue=-10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Minimum relative beta rotation (deg)") - # self.setPropertySettings("MinBetaRotation", condition) - # self.declareProperty(name="MaxBetaRotation", defaultValue=10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Maximum relative beta rotation (deg)") - # self.setPropertySettings("MaxBetaRotation", condition) - # - # # gamma rotation - # self.declareProperty(name="GammaRotation", defaultValue=False, - # doc="Refine rotation around third axis, gamma, for the components") - # condition = EnabledWhenProperty("GammaRotation", PropertyCriterion.IsNotDefault) - # self.declareProperty(name="MinGammaRotation", defaultValue=-10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Minimum relative gamma rotation (deg)") - # self.setPropertySettings("MinGammaRotation", condition) - # self.declareProperty(name="MaxGammaRotation", defaultValue=10.0, - # validator=FloatBoundedValidator(-90, 90), - # doc="Maximum relative gamma rotation (deg)") - # self.setPropertySettings("MaxGammaRotation", condition) - # - # properties = ["EulerConvention", "AlphaRotation", "MinAlphaRotation", "MaxAlphaRotation", - # "BetaRotation", "MinBetaRotation", "MaxBetaRotation", - # "GammaRotation", "MinGammaRotation", "MaxGammaRotation"] - # [self.setPropertyGroup(name, "Rotation") for name in properties] - # - # def validateInputs(self): - # """ - # Does basic validation for inputs - # """ - # issues = dict() - # - # peak_positions = self.getProperty('PeakPositions').value - # - # table_tof: TableWorkspace = self.getProperty('PeakCentersTofTable').value - # - # if 'detid' not in table_tof.getColumnNames(): - # issues['PeakCentersTofTable'] = 'PeakCentersTofTable is missing column "detid"' - # - # # The titles for table columns storing the TOF peak-center positions start with '@' - # column_names = [name for name in table_tof.getColumnNames() if name[0] == '@'] - # peak_count = len(peak_positions) # number of reference peak-center values - # if len(column_names) != peak_count: - # error_message = f'The number of table columns containing the peak center positions' \ - # f' {len(column_names)} is different than the number of peak positions {peak_count}' - # issues['PeakCentersTofTable'] = error_message - # - # # The titles for table columns storing the TOF peak-center positions do contain the values - # # of the reference peak-center positions in d-spacing units, up to a precision of 5 - # def with_precision(the_number, precision): - # r"""Analog of C++'s std::setprecision""" - # return round(the_number, precision - len(str(int(the_number)))) - # - # for column_name, peak_position in zip(column_names, sorted(peak_positions)): - # if (float(column_name[1:]) - with_precision(peak_position, 5)) > 1.e-5: - # issues['PeakCentersTofTable'] = f'{column_name} and {peak_position} differ up to precision 5' - # - # maskWS: MaskWorkspace = self.getProperty("MaskWorkspace").value - # if maskWS is not None: - # if maskWS.id() != 'MaskWorkspace': - # issues['MaskWorkspace'] = "MaskWorkspace must be empty or of type \"MaskWorkspace\"" - # # The mask workspace should contain as many spectra as rows in the TOFS table - # if maskWS.getNumberHistograms() != table_tof.rowCount(): - # error_message = 'The mask workspace must contain as many spectra as rows in the TOFS table' - # issues['MaskWorkspace'] = error_message - # - # # Need to get instrument in order to check if components are valid - # input_workspace = self.getProperty("InputWorkspace").value - # if bool(input_workspace) is True: - # wks_name = input_workspace.name() - # else: - # inputFilename = self.getProperty("InstrumentFilename").value - # if inputFilename == "": - # issues["InputWorkspace"] = "A Workspace or InstrumentFilename must be defined" - # return issues - # else: - # wks_name = "__alignedWorkspace" # a temporary workspace - # api.LoadEmptyInstrument(Filename=inputFilename, OutputWorkspace=wks_name) - # - # # Check if each component listed is defined in the instrument - # components = self.getProperty("ComponentList").value - # source_or_sample = self.getProperty("FitSourcePosition").value or self.getProperty("FitSamplePosition").value - # if len(components) <= 0 and not source_or_sample: - # issues['ComponentList'] = "Must supply components" - # else: - # get_component = api.mtd[wks_name].getInstrument().getComponentByName - # components = [component for component in components if get_component(component) is None] - # if len(components) > 0: - # issues['ComponentList'] = "Instrument has no component \"" \ - # + ','.join(components) + "\"" - # if wks_name == '__alignedWorkspace': - # api.DeleteWorkspace('__alignedWorkspace') # delete temporary workspace - # - # # This checks that something will actually be refined, - # if not (self.getProperty("Xposition").value - # or self.getProperty("Yposition").value - # or self.getProperty("Zposition").value - # or self.getProperty("AlphaRotation").value - # or self.getProperty("BetaRotation").value - # or self.getProperty("GammaRotation").value): - # issues["Xposition"] = "You must calibrate at least one position or rotation parameter." - # - # # Check that a position refinement is selected for sample/source - # if ((self.getProperty("FitSourcePosition").value - # or self.getProperty("FitSamplePosition").value) - # and not (self.getProperty("Xposition").value - # or self.getProperty("Yposition").value - # or self.getProperty("Zposition").value)): - # issues["Xposition"] = "If fitting source or sample, you must calibrate at least one position parameter." - # - # return issues - # - # # flake8: noqa: C901 - # def PyExec(self): - # table_tof = self.getProperty('PeakCentersTofTable').value - # self.peaks_tof = self._extract_tofs(table_tof) - # detector_count, peak_count = self.peaks_tof.shape - # table_tof = api.SortTableWorkspace(table_tof, Columns='detid') - # detID = table_tof.column('detid') - # peaks_ref = np.sort(self.getProperty('PeakPositions').value) # sort by increasing value - # self.peaks_ref = peaks_ref[np.newaxis, :] # shape = (1, peak_count) - # - # # Process input mask - # maskWS = self.getProperty("MaskWorkspace").value - # if maskWS is not None: - # mask = maskWS.extractY().flatten() # shape=(detector_count,) - # peaks_mask = np.tile(mask[:, np.newaxis], peak_count) # shape=(detector_count, peak_count) - # else: - # peaks_mask = np.zeros((detector_count, peak_count)) # no detectors are masked - # peaks_mask[np.isnan(self.peaks_tof)] = True - # # mask the defective detectors and missing peaks - # self.peaks_tof = np.ma.masked_array(self.peaks_tof, peaks_mask) - # - # input_workspace = self.getProperty('InputWorkspace').value - # - # adjustments_table_name = self.getProperty('AdjustmentsTable').value - # if len(adjustments_table_name) > 0: - # adjustments_table = self._initialize_adjustments_table(adjustments_table_name) - # saving_adjustments = True - # else: - # saving_adjustments = False - # - # self._eulerConvention = self.getProperty('EulerConvention').value - # - # output_workspace = self.getPropertyValue("OutputWorkspace") - # wks_name = '__alignedworkspace' # workspace whose counts will be DIFC values - # if bool(input_workspace) is True: - # api.CloneWorkspace(InputWorkspace=input_workspace, OutputWorkspace=wks_name) - # if output_workspace != str(input_workspace): - # api.CloneWorkspace(InputWorkspace=input_workspace, OutputWorkspace=output_workspace) - # else: - # api.LoadEmptyInstrument(Filename=self.getProperty("InstrumentFilename").value, OutputWorkspace=wks_name) - # - # # Make a dictionary of what options are being refined for sample/source. No rotation. - # for translation_option in self._optionsList[:3]: - # self._optionsDict[translation_option] = self.getProperty(translation_option).value - # for rotation_option in self._optionsList[3:]: - # self._optionsDict[rotation_option] = False - # - # # First fit L1 if selected for Source and/or Sample - # for component in "Source", "Sample": # fit first the source position, then the sample position - # if self.getProperty("Fit"+component+"Position").value: - # self._move = True - # if component == "Sample": - # comp = api.mtd[wks_name].getInstrument().getSample() - # else: - # comp = api.mtd[wks_name].getInstrument().getSource() - # componentName = comp.getFullName() - # logger.notice("Working on " + componentName + " Starting position is " + str(comp.getPos())) - # firstIndex = 0 - # lastIndex = detector_count - 1 - # - # self._initialPos = [comp.getPos().getX(), comp.getPos().getY(), comp.getPos().getZ(), - # 0, 0, 0] # no rotation - # - # # Set up x0 and bounds lists - # x0List = [] # initial X, Y, Z coordinates - # boundsList = [] # [(minX, maxX), (minZ, maxZ), (minZ, maxZ)] - # for iopt, translation_option in enumerate(self._optionsList[:3]): # iterate over X, Y, and Z - # if self._optionsDict[translation_option]: - # x0List.append(self._initialPos[iopt]) - # # default range for X is (x0 - 0.1m, x0 + 0.1m), same for Y and Z - # boundsList.append((self._initialPos[iopt] + self.getProperty("Min"+translation_option).value, - # self._initialPos[iopt] + self.getProperty("Max"+translation_option).value)) - # - # # scipy.opimize.minimize with the L-BFGS-B algorithm - # results: OptimizeResult = minimize(self._minimisation_func, x0=x0List, - # method='L-BFGS-B', - # args=(wks_name, componentName, firstIndex, lastIndex), - # bounds=boundsList) - # - # # Apply the results to the output workspace - # xmap = self._mapOptions(results.x) - # - # # Save translation and rotations, if requested - # if saving_adjustments: - # instrument = api.mtd[wks_name].getInstrument() - # name_finder = {'Source': instrument.getSource().getName(), - # 'Sample': instrument.getSample().getName()} - # component_adjustments = [name_finder[component]] + xmap[:3] + [0.0] * 4 # no rotations - # adjustments_table.addRow(component_adjustments) - # - # # Need to grab the component again, as things have changed - # kwargs = dict(X=xmap[0], Y=xmap[1], Z=xmap[2], RelativePosition=False, EnableLogging=False) - # api.MoveInstrumentComponent(wks_name, componentName, **kwargs) # adjust workspace - # api.MoveInstrumentComponent(output_workspace, componentName, **kwargs) # adjust workspace - # comp = api.mtd[wks_name].getInstrument().getComponentByName(componentName) - # logger.notice("Finished " + componentName + " Final position is " + str(comp.getPos())) - # self._move = False - # - # # Now fit all the remaining components, if any - # components = self.getProperty("ComponentList").value - # - # # Make a dictionary of what translational and rotational options are being refined. - # for opt in self._optionsList: - # self._optionsDict[opt] = self.getProperty(opt).value - # - # self._move = any([self._optionsDict[t] for t in ('Xposition', 'Yposition', 'Zposition')]) - # self._rotate = any([self._optionsDict[r] for r in ('AlphaRotation', 'BetaRotation', 'GammaRotation')]) - # - # prog = Progress(self, start=0, end=1, nreports=len(components)) - # get_component = api.mtd[wks_name].getInstrument().getComponentByName # shortcut - # for component in components: - # comp = get_component(component) - # firstDetID = self._getFirstDetID(comp) - # firstIndex = detID.index(firstDetID) # a row index in the input TOFS table - # lastDetID = self._getLastDetID(comp) - # lastIndex = detID.index(lastDetID) # a row index in the input TOFS table - # if lastDetID - firstDetID != lastIndex - firstIndex: - # raise RuntimeError("TOFS detid doesn't match instrument") - # - # eulerAngles: List[float] = comp.getRotation().getEulerAngles(self._eulerConvention) - # - # logger.notice("Working on " + comp.getFullName() + " Starting position is " + str(comp.getPos()) - # + " Starting rotation is " + str(eulerAngles)) - # - # x0List = [] - # self._initialPos = [comp.getPos().getX(), comp.getPos().getY(), comp.getPos().getZ(), - # eulerAngles[0], eulerAngles[1], eulerAngles[2]] - # - # boundsList = [] - # - # if np.all(peaks_mask[firstIndex:lastIndex + 1].astype(bool)): - # self.log().warning("All pixels in '%s' are masked. Skipping calibration." % component) - # continue - # - # for iopt, opt in enumerate(self._optionsList): - # if self._optionsDict[opt]: - # x0List.append(self._initialPos[iopt]) - # boundsList.append((self._initialPos[iopt] + self.getProperty("Min"+opt).value, - # self._initialPos[iopt] + self.getProperty("Max"+opt).value)) - # - # # scipy.opimize.minimize with the L-BFGS-B algorithm - # results: OptimizeResult = minimize(self._minimisation_func, x0=x0List, - # method='L-BFGS-B', - # args=(wks_name, component, firstIndex, lastIndex), - # bounds=boundsList) - # - # # Apply the results to the output workspace - # xmap = self._mapOptions(results.x) - # - # component_adjustments = [0.] * 7 # 3 for translation, 3 for rotation axis, 1 for rotation angle - # - # if self._move: - # kwargs = dict(X=xmap[0], Y=xmap[1], Z=xmap[2], RelativePosition=False, EnableLogging=False) - # api.MoveInstrumentComponent(wks_name, component, **kwargs) # adjust workspace - # api.MoveInstrumentComponent(output_workspace, component, **kwargs) # adjust workspace - # component_adjustments[:3] = xmap[:3] - # - # if self._rotate: - # (rotw, rotx, roty, rotz) = self._eulerToAngleAxis(xmap[3], xmap[4], xmap[5], self._eulerConvention) - # kwargs = dict(X=rotx, Y=roty, Z=rotz, Angle=rotw, RelativeRotation=False, EnableLogging=False) - # api.RotateInstrumentComponent(wks_name, component, **kwargs) # adjust workspace - # api.RotateInstrumentComponent(output_workspace, component, **kwargs) # adjust workspace - # component_adjustments[3:] = [rotx, roty, rotz, rotw] - # - # if saving_adjustments and (self._move or self._rotate): - # adjustments_table.addRow([component] + component_adjustments) - # - # # Need to grab the component object again, as things have changed - # comp = get_component(component) - # logger.notice("Finished " + comp.getFullName() + " Final position is " + str(comp.getPos()) - # + " Final rotation is " + str(comp.getRotation().getEulerAngles(self._eulerConvention))) - # - # prog.report() - # api.DeleteWorkspace(wks_name) - # self.setProperty("OutputWorkspace", output_workspace) - # logger.notice("Results applied to workspace "+wks_name) - # - # def _initialize_adjustments_table(self, table_name): - # r"""Create a table with appropriate column names for saving the adjustments to each component""" - # table = api.CreateEmptyTableWorkspace(OutputWorkspace=table_name) - # item_types = ['str', - # 'double', 'double', 'double', - # 'double', 'double', 'double', 'double'] - # for column_name, column_type in zip(self.adjustment_items, item_types): - # table.addColumn(name=column_name, type=column_type) - # return table - # - # def _extract_tofs(self, table_tofs: TableWorkspace) -> np.ndarray: - # r""" - # Extract the columns of the input table containing the peak centers, sorted by increasing value - # of the peak center in d-spacing units - # - # :param table_tofs: table of peak centers, in TOF units - # :return array of shape (detector_count, peak_count) - # """ - # # the title for the columns containing the peak centers begin with '@' - # indexes_and_titles = [(index, title) for index, title in enumerate(table_tofs.getColumnNames()) if '@' in title] - # column_indexes, titles = list(zip(*indexes_and_titles)) - # peak_tofs = np.array([table_tofs.column(i) for i in column_indexes]) # shape = (peak_count, detector_count) - # # sort by increasing peak centers (d-spacing units) - # peak_centers = np.array([float(title.replace('@', '')) for title in titles]) - # permutation = np.argsort(peak_centers) # reorder of indices guarantee increase in d-spacing - # peak_tofs = peak_tofs[permutation] # sort by increasing d-spacing - # return np.transpose(peak_tofs) # shape = (detector_count, peak_count) - # - # def _minimisation_func(self, x_0, wks_name, component, firstIndex, lastIndex): - # """ - # Basic minimization function used. Returns the sum of the absolute values for the fractional peak - # deviations: - # - # .. math:: - # - # \\sum_i^{N_d}\\sum_j^{N_p} (1 - m_{i,j}) \\frac{|d_{i,j} - d_j^*|}{d_j^*} - # - # where :math:`N_d` is the number of detectors in the bank, :math:`N_p` is the number of reference peaks, and - # :math:`m_{i,j}` is the mask for peak :math:`j` and detector :math:`i`. The mask evaluates to 1 if the - # detector is defective or the peak is missing in the detector, otherwise the mask evaluates to zero. - # - # There's an implicit one-to-correspondence between array index of ``difc`` and workspace index of ``wks_name``, - # that is, between row index of the input TOFS table and workspace index of ``wks_name``. - # - # @param x_0 :: list of length 3 (new XYZ coordinates of the component) or length 6 (XYZ and rotation coords) - # @param wks_name :: name of a workspace with an embedded instrument. The instrument will be adjusted according to - # the new coordinates ``x_0`` for instrument component ``component``. It's pixel spectra will contain the new DIFC - # @param component :: name of the instrument component to be optimized - # @param firstIndex :: workspace index of first index of ``difc`` array to be considered when comparing old - # and new DIFC values. When fitting the source or sample, this is the first spectrum index. - # @param lastIndex :: workspace index of last index of ``difc`` array to be considered when comparing old - # and new DIFC values. When fitting the source or sample, this is the last row number of the input - # TOFS table. - # - # @return Chi-square value between old and new DIFC values for the unmasked spectra - # """ - # xmap = self._mapOptions(x_0) # pad null rotations when x_0 contains only translations - # - # if self._move: - # api.MoveInstrumentComponent(wks_name, component, X=xmap[0], Y=xmap[1], Z=xmap[2], - # RelativePosition=False, EnableLogging=False) - # - # if self._rotate: - # (rotw, rotx, roty, rotz) = self._eulerToAngleAxis(xmap[3], xmap[4], xmap[5], self._eulerConvention) # YZX - # api.RotateInstrumentComponent(wks_name, component, X=rotx, Y=roty, Z=rotz, Angle=rotw, - # RelativeRotation=False, EnableLogging=False) - # - # api.CalculateDIFC(InputWorkspace=wks_name, OutputWorkspace=wks_name, EnableLogging=False) - # difc = api.mtd[wks_name].extractY().flatten()[firstIndex: lastIndex + 1] - # peaks_d = self.peaks_tof[firstIndex: lastIndex + 1] / difc[:, np.newaxis] # peak centers in d-spacing units - # - # # calculate the fractional peak center deviations, then sum their absolute values - # return np.sum(np.abs((peaks_d - self.peaks_ref) / self.peaks_ref)) - # - # def _getFirstDetID(self, component): - # """ - # recursive search to find first detID of a component - # - # @param component :: reference to a detector component object - # - # @returns detector ID (``int``) of the first detector in the component - # """ - # if component.type() == 'DetectorComponent' or component.type() == 'GridDetectorPixel': - # return component.getID() - # else: - # return self._getFirstDetID(component[0]) - # - # def _getLastDetID(self, component): - # """ - # recursive search to find last detID of a component - # - # @param component :: reference to a detector component object - # - # @returns detector ID (`int`) of the last detector in the component """ - # if component.type() == 'DetectorComponent' or component.type() == 'GridDetectorPixel': - # return component.getID() - # else: - # return self._getLastDetID(component[component.nelements() - 1]) - # - # def _mapOptions(self, inX): - # """ - # Creates an array combining the refining and constant variables - # This is required because scipy.optimise.minimize expects a constant - # number of variables, so need to be able to maps any number of - # inputs to six outputs. - # - # @param inX :: list of length 3 or 6 - # - # @return list of length 6 - # """ - # x0_index = 0 - # out = [] - # for opt in self._optionsList: - # if self._optionsDict[opt]: - # out.append(inX[x0_index]) - # x0_index += 1 - # else: - # out.append(self._initialPos[self._optionsList.index(opt)]) - # return out - # - # def _eulerToQuat(self, alpha, beta, gamma, convention): - # """ - # Convert Euler angles to a quaternion - # """ - # getV3D = {'X': V3D(1, 0, 0), 'Y': V3D(0, 1, 0), 'Z': V3D(0, 0, 1)} - # return (Quat(alpha, getV3D[convention[0]]) * Quat(beta, getV3D[convention[1]]) - # * Quat(gamma, getV3D[convention[2]])) - # - # def _eulerToAngleAxis(self, alpha, beta, gamma, convention): - # """ - # Find the Euler axis and Euler angle - # - # @param alpha :: rotation angle, in degrees, around the first axis of `convention` - # @param beta :: rotation angle, in degrees, around the second axis of `convention` - # @param gamma :: rotation angle, in degrees, around the third axis of `convention` - # @param convention :: string, e.g. 'YZX' - # - # @return Euler angle, and three direct cosines defining the Euler axis - # """ - # quat = self._eulerToQuat(alpha, beta, gamma, convention) - # if quat[0] == 1: - # return 0, 0, 0, 1 - # deg = math.acos(quat[0]) - # scale = math.sin(deg) - # deg *= 360.0 / math.pi - # ax0 = quat[1] / scale - # ax1 = quat[2] / scale - # ax2 = quat[3] / scale - # return deg, ax0, ax1, ax2 - # - # - # try: - # from scipy.optimize import minimize, OptimizeResult - # AlgorithmFactory.subscribe(AlignComponents) - # except ImportError: - # logger.debug('Failed to subscribe algorithm AlignComponets; cannot import minimize from scipy.optimize')eak centers (d-spacing units) peak_centers = np.array([float(title.replace('@', '')) for title in titles]) permutation = np.argsort(peak_centers) # reorder of indices guarantee increase in d-spacing peak_tofs = peak_tofs[permutation] # sort by increasing d-spacing