From e1dd11881c2c0c03ec5188462af703c7fe9e2019 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 26 Apr 2023 12:14:23 +0100 Subject: [PATCH 01/11] add functionality for different fit type options --- .../calibrate_vesuvio_analysis.py | 86 ++++++++----- .../calibrate_vesuvio_fit.py | 119 +++++++++++------- 2 files changed, 129 insertions(+), 76 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py index f439f60f..9d64c8af 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py @@ -137,7 +137,7 @@ def PyExec(self): self._invalid_detectors.add_invalid_detectors(invalid_detectors) E1_peak_fits_back = mtd[self._current_workspace + '_E1_back_Peak_Parameters'].getNames() - self._calculate_final_energy(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type != "Individual") + self._calculate_final_energy(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type) # calibrate L1 for backscattering detectors based on the averaged E1 value and calibrated theta values self._calculate_final_flight_path(E1_peak_fits_back[0], EVSGlobals.BACKSCATTERING_RANGE) @@ -315,7 +315,7 @@ def _calculate_scattering_angle(self, table_name, spec_list): self._set_table_column(self._current_workspace, 'theta', masked_theta, spec_list) self._set_table_column(self._current_workspace, 'theta_Err', theta_error, spec_list) - def _calculate_final_energy(self, peak_table, spec_list, calculate_global): + def _calculate_final_energy(self, peak_table, spec_list, fit_type): """ Calculate the final energy using the fitted peak centres of a run. @@ -324,12 +324,17 @@ def _calculate_final_energy(self, peak_table, spec_list, calculate_global): """ spec_range = EVSGlobals.DETECTOR_RANGE[1] + 1 - EVSGlobals.DETECTOR_RANGE[0] - mean_E1 = np.empty(spec_range) - E1_error = np.empty(spec_range) - global_E1 = np.empty(spec_range) - global_E1_error = np.empty(spec_range) if not self._E1_value_and_error: + + if fit_type == 'Both': + indv_peak_table = peak_table[0] + glob_peak_table = peak_table[1] + elif fit_type == 'Shared': + glob_peak_table = peak_table[0] + else: + indv_peak_table = peak_table[0] + t0 = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list) L0 = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list) @@ -337,44 +342,59 @@ def _calculate_final_energy(self, peak_table, spec_list, calculate_global): theta = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list) r_theta = EVSMiscFunctions.calculate_r_theta(self._sample_mass, theta) - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, peak_table[0]) + if fit_type != 'Shared': + mean_E1 = np.empty(spec_range) + E1_error = np.empty(spec_range) - delta_t = (peak_centres - t0) / 1e+6 - v1 = (L0 * r_theta + L1) / delta_t + peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) - E1 = 0.5 * scipy.constants.m_n * v1 ** 2 - E1 /= EVSGlobals.MEV_CONVERSION + delta_t = (peak_centres - t0) / 1e+6 + v1 = (L0 * r_theta + L1) / delta_t - mean_E1_val = np.nanmean(E1) - E1_error_val = np.nanstd(E1) + E1 = 0.5*scipy.constants.m_n*v1**2 + E1 /= EVSGlobals.MEV_CONVERSION - else: - mean_E1_val = self._E1_value_and_error[0] - E1_error_val = self._E1_value_and_error[1] + mean_E1_val = np.nanmean(E1) + E1_error_val = np.nanstd(E1) + + mean_E1.fill(mean_E1_val) + E1_error.fill(E1_error_val) + + self._set_table_column(self._current_workspace, 'E1', mean_E1) + self._set_table_column(self._current_workspace, 'E1_Err', E1_error) + + if fit_type != 'Individual': + global_E1 = np.empty(spec_range) + global_E1_error = np.empty(spec_range) - mean_E1.fill(mean_E1_val) - E1_error.fill(E1_error_val) + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(glob_peak_table, 'f1.LorentzPos', [0]) + peak_centre = [peak_centre] * spec_range - self._set_table_column(self._current_workspace, 'E1', mean_E1) - self._set_table_column(self._current_workspace, 'E1_Err', E1_error) + delta_t = (peak_centre - t0) / 1e+6 + v1 = (L0 * r_theta + L1) / delta_t + E1 = 0.5*scipy.constants.m_n*v1**2 + E1 /= EVSGlobals.MEV_CONVERSION - if calculate_global: # This fn will need updating for the only global option - peak_centre = EVSMiscFunctions.read_fitting_result_table_column(peak_table[1], 'f1.LorentzPos', [0]) - peak_centre = [peak_centre] * len(peak_centres) + global_E1_val = np.nanmean(E1) + global_E1_error_val = np.nanstd(E1) - delta_t = (peak_centre - t0) / 1e+6 - v1 = (L0 * r_theta + L1) / delta_t - E1 = 0.5 * scipy.constants.m_n * v1 ** 2 - E1 /= EVSGlobals.MEV_CONVERSION + global_E1.fill(global_E1_val) + global_E1_error.fill(global_E1_error_val) - global_E1_val = np.nanmean(E1) - global_E1_error_val = np.nanstd(E1) + self._set_table_column(self._current_workspace, 'global_E1', global_E1) + self._set_table_column(self._current_workspace, 'global_E1_Err', global_E1_error) + else: + mean_E1 = np.empty(spec_range) + E1_error = np.empty(spec_range) + + mean_E1_val = self._E1_value_and_error[0] + E1_error_val = self._E1_value_and_error[1] - global_E1.fill(global_E1_val) - global_E1_error.fill(global_E1_error_val) + mean_E1.fill(mean_E1_val) + E1_error.fill(E1_error_val) - self._set_table_column(self._current_workspace, 'global_E1', global_E1) - self._set_table_column(self._current_workspace, 'global_E1_Err', global_E1_error) + self._set_table_column(self._current_workspace, 'E1', mean_E1) + self._set_table_column(self._current_workspace, 'E1_Err', E1_error) def _setup(self): """ diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index efef766a..35ea9084 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -555,60 +555,47 @@ def _fit_peaks(self): """ estimated_peak_positions_all_peaks = self._estimate_peak_positions() - num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape + if self._shared_parameter_fit_type != "Shared": + num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape - self._prog_reporter = Progress(self, 0.0, 1.0, num_estimated_peaks * num_spectra) + self._prog_reporter = Progress(self, 0.0, 1.0, num_estimated_peaks*num_spectra) - self._output_parameter_tables = [] - self._peak_fit_workspaces = [] + self._output_parameter_tables = [] + self._peak_fit_workspaces = [] + for peak_index, estimated_peak_positions in enumerate(estimated_peak_positions_all_peaks): - for peak_index, estimated_peak_positions in enumerate(estimated_peak_positions_all_peaks): + self._peak_fit_workspaces_by_spec = [] + output_parameter_table_name = self._output_workspace_name + '_Peak_%d_Parameters' % peak_index + output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name) + for spec_index, peak_position in enumerate(estimated_peak_positions): + fit_workspace_name = self._fit_peak(peak_index, spec_index, peak_position, output_parameter_table_name, + output_parameter_table_headers) + self._peak_fit_workspaces_by_spec.append(fit_workspace_name) - self._peak_fit_workspaces_by_spec = [] - output_parameter_table_name = self._output_workspace_name + '_Peak_%d_Parameters' % peak_index - output_parameter_table_headers = self._create_parameter_table_and_output_headers( - output_parameter_table_name) - for spec_index, peak_position in enumerate(estimated_peak_positions): - fit_workspace_name = self._fit_peak(peak_index, spec_index, peak_position, output_parameter_table_name, - output_parameter_table_headers) - self._peak_fit_workspaces_by_spec.append(fit_workspace_name) + self._output_parameter_tables.append(output_parameter_table_name) + self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) - self._output_parameter_tables.append(output_parameter_table_name) - self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) - - GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') + GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') if self._shared_parameter_fit_type != "Individual": - self._shared_parameter_fit(output_parameter_table_name, output_parameter_table_headers) - - def _shared_parameter_fit(self, output_parameter_table_name, output_parameter_table_headers): - init_Gaussian_FWHM = EVSMiscFunctions.read_fitting_result_table_column(output_parameter_table_name, 'f1.GaussianFWHM', - self._spec_list) - init_Gaussian_FWHM = np.nanmean(init_Gaussian_FWHM[init_Gaussian_FWHM != 0]) - init_Lorentz_FWHM = EVSMiscFunctions.read_fitting_result_table_column(output_parameter_table_name, 'f1.LorentzFWHM', - self._spec_list) - init_Lorentz_FWHM = np.nanmean(init_Lorentz_FWHM[init_Lorentz_FWHM != 0]) - init_Lorentz_Amp = EVSMiscFunctions.read_fitting_result_table_column(output_parameter_table_name, 'f1.LorentzAmp', - self._spec_list) - init_Lorentz_Amp = np.nanmean(init_Lorentz_Amp[init_Lorentz_Amp != 0]) - init_Lorentz_Pos = EVSMiscFunctions.read_fitting_result_table_column(output_parameter_table_name, 'f1.LorentzPos', - self._spec_list) - init_Lorentz_Pos = np.nanmean(init_Lorentz_Pos[init_Lorentz_Pos != 0]) - - initial_params = {'A0': 0.0, 'A1': 0.0, 'LorentzAmp': init_Lorentz_Amp, 'LorentzPos': init_Lorentz_Pos, - 'LorentzFWHM': init_Lorentz_FWHM, 'GaussianFWHM': init_Gaussian_FWHM} - - self._invalid_detectors.identify_and_set_invalid_detectors_from_range(self._spec_list, output_parameter_table_name) - - self._fit_shared_parameter(self._invalid_detectors.get_invalid_detectors_index(self._spec_list), initial_params, - output_parameter_table_headers) + estimated_peak_position = np.mean(estimated_peak_positions_all_peaks) + output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_Parameters' + output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name) + + self._fit_shared_peak(self._spec_list[0], estimated_peak_position, output_parameter_table_name, + output_parameter_table_headers) + + if self._shared_parameter_fit_type == 'Both': + mtd[self._output_workspace_name + '_Peak_Parameters'].add(output_parameter_table_name) + else: + GroupWorkspaces(output_parameter_table_name, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_table_name, output_parameter_table_headers): spec_number = self._spec_list[0] + spec_index self._prog_reporter.report("Fitting peak %d to spectrum %d" % (peak_index, spec_number)) - peak_params = self._find_peaks_and_output_params(peak_index, spec_number, spec_index, peak_position) + peak_params = self._find_peaks_and_output_params(spec_number, peak_position, peak_index, spec_index) fit_func_string = self._build_function_string(peak_params) xmin, xmax = self._find_fit_x_window(peak_params) fit_output_name = '__' + self._output_workspace_name + '_Peak_%d_Spec_%d' % (peak_index, spec_index) @@ -624,9 +611,55 @@ def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_tabl self._del_fit_workspaces(ncm, fit_params, fws) return fit_workspace_name + + def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers): + peak_params = self._find_peaks_and_output_params(spec_range, peak_position) + fit_func = self._build_function_string(peak_params) + start_str = 'composite=MultiDomainFunction, NumDeriv=1;' + validSpecs = mtd[self._sample] + # need to decide how to remove invalid spectra from multifit + # MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) + # validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') + n_valid_specs = validSpecs.getNumberHistograms() + + fit_func = ('(composite=CompositeFunction, NumDeriv=false, $domains=i;' + fit_func[:-1] + ');') * n_valid_specs + composite_func = start_str + fit_func[:-1] + + ties = ','.join(f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1,n_valid_specs)) + func_string = composite_func + f';ties=({ties})' + xmin, xmax = self._find_fit_x_window(peak_params) + fit_output_name = '__' + self._output_workspace_name + '_Peak_0' - def _find_peaks_and_output_params(self, peak_index, spec_number, spec_index, peak_position): - peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index) + # create new workspace for each spectra + x = validSpecs.readX(0) + y = validSpecs.readY(0) + e = validSpecs.readE(0) + out_ws = self._sample + '_Spec_0' + CreateWorkspace(DataX=x, DataY=y, DataE=e, NSpec=1, OutputWorkspace=out_ws) + + other_inputs = [CreateWorkspace(DataX=validSpecs.readX(j), DataY=validSpecs.readY(j), DataE=validSpecs.readE(j), + OutputWorkspace=f'{self._sample}_Spec_{j}') + for j in range(1,n_valid_specs)] + + added_args = {f'InputWorkspace_{i + 1}': v for i,v in enumerate(other_inputs)} + + status, chi2, ncm, fit_params, fws, func, cost_func = Fit(Function=func_string, InputWorkspace=out_ws, IgnoreInvalidData=True, + StartX=xmin, EndX=xmax, + CalcErrors=True, Output=fit_output_name, + Minimizer='SteepestDescent,RelError=1e-8', **added_args) + [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0,n_valid_specs)] + + output_headers = ['f0.'+ name for name in output_parameter_table_headers] + + self._output_fit_params_to_table_ws(0, fit_params, output_parameter_table_name, + output_headers) + self._del_fit_workspaces(ncm, fit_params, fws) + + def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index=None, spec_index=None): + if spec_index and peak_index: + peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index) + else: + peak_table_name = '__' + self._sample + '_peak_table' find_peak_params = self._get_find_peak_parameters(spec_number, [peak_position]) FindPeaks(InputWorkspace=self._sample, WorkspaceIndex=spec_index, PeaksList=peak_table_name, **find_peak_params) if mtd[peak_table_name].rowCount() == 0: From e877ba7fa1d2505b273e30f4fbf8bd1e22ed794f Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 26 Apr 2023 14:03:24 +0100 Subject: [PATCH 02/11] remove function as no longer needed --- .../calibrate_vesuvio_fit.py | 73 ------------------- 1 file changed, 73 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index 35ea9084..c4082e74 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -696,79 +696,6 @@ def _del_fit_workspaces(self, ncm, fit_params, fws): if not self._create_output: DeleteWorkspace(fws) - def _fit_shared_parameter(self, invalid_spectra, initial_params, param_names): - """ - Fit peaks to all detectors, with one set of fit parameters for all detectors. - """ - - shared_peak_table = self._output_workspace_name + '_shared_peak_parameters' - CreateEmptyTableWorkspace(OutputWorkspace=shared_peak_table) - - err_names = [name + '_Err' for name in param_names] - col_names = [element for tupl in zip(param_names, err_names) for element in tupl] - - mtd[shared_peak_table].addColumn('int', 'Spectrum') - for name in col_names: - mtd[shared_peak_table].addColumn('double', name) - - fit_func = self._build_function_string(initial_params) - start_str = 'composite=MultiDomainFunction, NumDeriv=1;' - sample_ws = mtd[self._sample] - MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) - validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') - n_valid_specs = validSpecs.getNumberHistograms() - - fit_func = ('(composite=CompositeFunction, NumDeriv=false, $domains=i;' + fit_func[:-1] + ');') * n_valid_specs - composite_func = start_str + fit_func[:-1] - - ties = ','.join( - f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1, n_valid_specs)) - func_string = composite_func + f';ties=({ties})' - - fit_output_name = '__' + self._output_workspace_name + '_Peak_0' - - xmin, xmax = None, None - - # create new workspace for each spectra - x = validSpecs.readX(0) - y = validSpecs.readY(0) - e = validSpecs.readE(0) - out_ws = self._sample + '_Spec_0' - CreateWorkspace(DataX=x, DataY=y, DataE=e, NSpec=1, OutputWorkspace=out_ws) - - other_inputs = [CreateWorkspace(DataX=validSpecs.readX(j), DataY=validSpecs.readY(j), DataE=validSpecs.readE(j), - OutputWorkspace=f'{self._sample}_Spec_{j}') - for j in range(1, n_valid_specs)] - - added_args = {f'InputWorkspace_{i + 1}': v for i, v in enumerate(other_inputs)} - - print('starting global fit') - - status, chi2, ncm, params, fws, func, cost_func = Fit(Function=func_string, InputWorkspace=out_ws, - IgnoreInvalidData=True, - StartX=xmin, EndX=xmax, - CalcErrors=True, Output=fit_output_name, - Minimizer='SteepestDescent,RelError=1e-8', **added_args) - [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0, n_valid_specs)] - - fit_values = dict(zip(params.column(0), params.column(1))) - fit_errors = dict(zip(params.column(0), params.column(2))) - - row_values = [fit_values['f0.' + name] for name in param_names] - row_errors = [fit_errors['f0.' + name] for name in param_names] - row = [element for tupl in zip(row_values, row_errors) for element in tupl] - - mtd[shared_peak_table].addRow([0] + row) - - DeleteWorkspace(ncm) - DeleteWorkspace(params) - if not self._create_output: - DeleteWorkspace(fws) - - DeleteWorkspace(validSpecs) - - mtd[self._output_workspace_name + '_Peak_Parameters'].add(shared_peak_table) - def _get_find_peak_parameters(self, spec_number, peak_centre, unconstrained=False): """ Get find peak parameters From a26051caf8af509e08acf7d46fa6eacba8a2073a Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 26 Apr 2023 15:14:06 +0100 Subject: [PATCH 03/11] fix unit test failures --- .../calibrate_vesuvio_fit.py | 12 +++--- .../tests/unit/test_calibrate_vesuvio_fit.py | 37 ++++++++----------- 2 files changed, 22 insertions(+), 27 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index c4082e74..f5ef748c 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -572,10 +572,11 @@ def _fit_peaks(self): output_parameter_table_headers) self._peak_fit_workspaces_by_spec.append(fit_workspace_name) - self._output_parameter_tables.append(output_parameter_table_name) - self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) + self._output_parameter_tables.append(output_parameter_table_name) + self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) - GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') + if self._shared_parameter_fit_type == 'Individual': + GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') if self._shared_parameter_fit_type != "Individual": estimated_peak_position = np.mean(estimated_peak_positions_all_peaks) @@ -586,7 +587,8 @@ def _fit_peaks(self): output_parameter_table_headers) if self._shared_parameter_fit_type == 'Both': - mtd[self._output_workspace_name + '_Peak_Parameters'].add(output_parameter_table_name) + self._output_parameter_tables.append(output_parameter_table_name) + GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') else: GroupWorkspaces(output_parameter_table_name, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') @@ -656,7 +658,7 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam self._del_fit_workspaces(ncm, fit_params, fws) def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index=None, spec_index=None): - if spec_index and peak_index: + if spec_index is not None and peak_index is not None: peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index) else: peak_table_name = '__' + self._sample + '_peak_table' diff --git a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py index 3c8bce65..04a7d572 100644 --- a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py @@ -660,39 +660,32 @@ def test_fit_peaks_individual(self, group_workspaces_mock): group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters') - @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._shared_parameter_fit') + @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak') @patch('calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces') - def test_fit_peaks_shared(self, group_workspaces_mock, shared_parameter_fit_mock): + def test_fit_peaks_shared(self, group_workspaces_mock, fit_shared_peak_mock): alg = EVSCalibrationFit() alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]])) alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') alg._output_workspace_name = 'output_ws_name' alg._shared_parameter_fit_type = 'Shared' + alg._spec_list = [3,4,5] alg._fit_peaks() - alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Peak_0_Parameters'), - call('output_ws_name_Peak_1_Parameters')]) - alg._fit_peak.assert_has_calls([call(0, 0, 5, 'output_ws_name_Peak_0_Parameters', ['a', 'b', 'c']), - call(0, 1, 10, 'output_ws_name_Peak_0_Parameters', ['a', 'b', 'c']), - call(0, 2, 15, 'output_ws_name_Peak_0_Parameters', ['a', 'b', 'c']), - call(1, 0, 2.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']), - call(1, 1, 7.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']), - call(1, 2, 10.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c'])]) - self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2']], - alg._peak_fit_workspaces) - group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters'], + alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Shared_Peak_Parameters')]) + group_workspaces_mock.assert_called_once_with('output_ws_name_Shared_Peak_Parameters', OutputWorkspace='output_ws_name_Peak_Parameters') - shared_parameter_fit_mock.assert_called_once_with('output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']) + fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c']) - @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._shared_parameter_fit') + @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak') @patch('calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces') - def test_fit_peaks_both(self, group_workspaces_mock, shared_parameter_fit_mock): + def test_fit_peaks_both(self, group_workspaces_mock, fit_shared_peak_mock): alg = EVSCalibrationFit() alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]])) alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') alg._output_workspace_name = 'output_ws_name' alg._shared_parameter_fit_type = 'Both' + alg._spec_list = [3,4,5] alg._fit_peaks() alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Peak_0_Parameters'), call('output_ws_name_Peak_1_Parameters')]) @@ -704,9 +697,9 @@ def test_fit_peaks_both(self, group_workspaces_mock, shared_parameter_fit_mock): call(1, 2, 10.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c'])]) self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2']], alg._peak_fit_workspaces) - group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters'], - OutputWorkspace='output_ws_name_Peak_Parameters') - shared_parameter_fit_mock.assert_called_once_with('output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']) + group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters', + 'output_ws_name_Shared_Peak_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters') + fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c']) @staticmethod def _setup_alg_mocks_fit_peak(): @@ -734,7 +727,7 @@ def test_fit_peak(self, mock_fit): fit_workspace_name = alg._fit_peak(0, 0, 5, 'output_Peak_0_Parameters', ['a', 'b', 'c']) - alg._find_peaks_and_output_params.assert_called_once_with(0, 30, 0, 5) + alg._find_peaks_and_output_params.assert_called_once_with(30, 5, 0, 0) alg._build_function_string.assert_called_once_with('peak_params') alg._find_fit_x_window.assert_called_once_with('peak_params') mock_fit.assert_called_once_with(Function='function_string', InputWorkspace='sample', @@ -759,7 +752,7 @@ def test_find_peaks_and_output_params(self, find_peaks_mock, delete_workspace_mo peak_table_name_ws.rowCount.return_value = 5 mtd_mock.__getitem__.return_value = peak_table_name_ws - alg._find_peaks_and_output_params(0, 30, 3, 5) + alg._find_peaks_and_output_params(30, 5, 0, 3) alg._get_find_peak_parameters.assert_called_once_with(30, [5]) find_peaks_mock.assert_called_once_with(InputWorkspace='sample', WorkspaceIndex=3, PeaksList='__sample_peaks_table_0_3', @@ -784,7 +777,7 @@ def test_find_peaks_and_output_params_no_peaks_found(self, find_peak_params_mock sys_mock.exit.side_effect = Exception("Emulate Sys Exit") with self.assertRaises(Exception): - alg._find_peaks_and_output_params(2, 32, 2, 7.5) + alg._find_peaks_and_output_params(32, 7.5, 2, 2) find_peak_params_mock.assert_called_once_with(32, [7.5]) find_peaks_mock.assert_called_once_with(InputWorkspace='sample', WorkspaceIndex=2, PeaksList='__sample_peaks_table_2_2', From 6c33460b84d7dd38c5160aca582cd49fff59f06a Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Thu, 27 Apr 2023 10:20:46 +0100 Subject: [PATCH 04/11] use lm to fit instead of steepest descent --- .../calibration_scripts/calibrate_vesuvio_fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index f5ef748c..d85bdbb7 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -648,7 +648,7 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam status, chi2, ncm, fit_params, fws, func, cost_func = Fit(Function=func_string, InputWorkspace=out_ws, IgnoreInvalidData=True, StartX=xmin, EndX=xmax, CalcErrors=True, Output=fit_output_name, - Minimizer='SteepestDescent,RelError=1e-8', **added_args) + Minimizer='Levenberg-Marquardt,RelError=1e-8', **added_args) [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0,n_valid_specs)] output_headers = ['f0.'+ name for name in output_parameter_table_headers] From 66d0f1ba40cde87b9bfac671f59742d3d9e44194 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Thu, 27 Apr 2023 15:53:03 +0000 Subject: [PATCH 05/11] allow fitting of multiple peaks and add new system test --- .../calibrate_vesuvio_analysis.py | 56 ++++++++++++++----- .../calibrate_vesuvio_fit.py | 40 +++++++------ .../tests/system/test_system_analysis.py | 26 ++++++++- .../tests/unit/test_calibrate_vesuvio_fit.py | 27 +++++---- 4 files changed, 104 insertions(+), 45 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py index 9d64c8af..f4259e56 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py @@ -140,7 +140,7 @@ def PyExec(self): self._calculate_final_energy(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type) # calibrate L1 for backscattering detectors based on the averaged E1 value and calibrated theta values - self._calculate_final_flight_path(E1_peak_fits_back[0], EVSGlobals.BACKSCATTERING_RANGE) + self._calculate_final_flight_path(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type) # calibrate L1 for frontscattering detectors based on the averaged E1 value and calibrated theta values E1_fit_front = self._current_workspace + '_E1_front' @@ -153,7 +153,7 @@ def PyExec(self): self._invalid_detectors.add_invalid_detectors(invalid_detectors) E1_peak_fits_front = mtd[self._current_workspace + '_E1_front_Peak_Parameters'].getNames() - self._calculate_final_flight_path(E1_peak_fits_front[0], EVSGlobals.FRONTSCATTERING_RANGE) + self._calculate_final_flight_path(E1_peak_fits_front, EVSGlobals.FRONTSCATTERING_RANGE, self._shared_parameter_fit_type) # make the fitted parameters for this iteration the input to the next iteration. table_group.append(self._current_workspace) @@ -245,33 +245,61 @@ def _calculate_incident_flight_path(self, table_name, spec_list): DeleteWorkspace(L0_param_table) - def _calculate_final_flight_path(self, peak_table, spec_list): + def _calculate_final_flight_path(self, peak_table, spec_list, fit_type): """ Calculate the final flight path using the values for energy. This also uses the old value for L1 loaded from the parameter file. @param spec_list - spectrum range to calculate t0 for. """ + if fit_type == 'Both': + indv_peak_table = peak_table[0] + glob_peak_table = peak_table[1] + E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list) + E1 *= EVSGlobals.MEV_CONVERSION + global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list) + global_E1 *= EVSGlobals.MEV_CONVERSION + elif fit_type == 'Shared': + glob_peak_table = peak_table[0] + global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list) + global_E1 *= EVSGlobals.MEV_CONVERSION + else: + indv_peak_table = peak_table[0] + E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list) + E1 *= EVSGlobals.MEV_CONVERSION - E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list) t0 = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list) t0_error = EVSMiscFunctions.read_table_column(self._current_workspace, 't0_Err', spec_list) L0 = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list) theta = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list) r_theta = EVSMiscFunctions.calculate_r_theta(self._sample_mass, theta) - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, peak_table) + if fit_type != 'Shared': + + peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) + + delta_t = (peak_centres - t0) / 1e+6 + delta_t_error = t0_error / 1e+6 + v1 = np.sqrt( 2 * E1 /scipy.constants.m_n) + L1 = v1 * delta_t - L0 * r_theta + L1_error = v1 * delta_t_error - delta_t = (peak_centres - t0) / 1e+6 - delta_t_error = t0_error / 1e+6 + self._set_table_column(self._current_workspace, 'L1', L1, spec_list) + self._set_table_column(self._current_workspace, 'L1_Err', L1_error, spec_list) - E1 *= EVSGlobals.MEV_CONVERSION - v1 = np.sqrt( 2 * E1 /scipy.constants.m_n) - L1 = v1 * delta_t - L0 * r_theta - L1_error = v1 * delta_t_error + if fit_type != 'Individual': + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(glob_peak_table, 'f1.LorentzPos', [0]) + peak_centre_list = np.empty(spec_list[1] + 1 -spec_list[0]) + peak_centre_list.fill(float(peak_centre)) - self._set_table_column(self._current_workspace, 'L1', L1, spec_list) - self._set_table_column(self._current_workspace, 'L1_Err', L1_error, spec_list) + delta_t = (peak_centre_list - t0) / 1e+6 + delta_t_error = t0_error / 1e+6 + global_v1 = np.sqrt( 2 * global_E1 /scipy.constants.m_n) + global_L1 = global_v1 * delta_t - L0 * r_theta + global_L1_error = global_v1 * delta_t_error + + self._set_table_column(self._current_workspace, 'global_L1', global_L1, spec_list) + self._set_table_column(self._current_workspace, 'global_L1_Err', global_L1_error, spec_list) def _calculate_scattering_angle(self, table_name, spec_list): """ @@ -383,7 +411,7 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): self._set_table_column(self._current_workspace, 'global_E1', global_E1) self._set_table_column(self._current_workspace, 'global_E1_Err', global_E1_error) - else: + else: mean_E1 = np.empty(spec_range) E1_error = np.empty(spec_range) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index d85bdbb7..bc534967 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -555,13 +555,12 @@ def _fit_peaks(self): """ estimated_peak_positions_all_peaks = self._estimate_peak_positions() + num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape + self._output_parameter_tables = [] + self._peak_fit_workspaces = [] if self._shared_parameter_fit_type != "Shared": - num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape - self._prog_reporter = Progress(self, 0.0, 1.0, num_estimated_peaks*num_spectra) - self._output_parameter_tables = [] - self._peak_fit_workspaces = [] for peak_index, estimated_peak_positions in enumerate(estimated_peak_positions_all_peaks): self._peak_fit_workspaces_by_spec = [] @@ -579,18 +578,21 @@ def _fit_peaks(self): GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') if self._shared_parameter_fit_type != "Individual": - estimated_peak_position = np.mean(estimated_peak_positions_all_peaks) - output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_Parameters' - output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name) + for peak_index, estimated_peak_position in enumerate(np.mean(estimated_peak_positions_all_peaks,1)): - self._fit_shared_peak(self._spec_list[0], estimated_peak_position, output_parameter_table_name, - output_parameter_table_headers) + output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_%d_Parameters' % peak_index + output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name) + + fit_workspace_name = self._fit_shared_peak(peak_index, self._spec_list[0], estimated_peak_position, output_parameter_table_name, + output_parameter_table_headers) - if self._shared_parameter_fit_type == 'Both': self._output_parameter_tables.append(output_parameter_table_name) + self._peak_fit_workspaces.append(fit_workspace_name) + + if self._shared_parameter_fit_type == 'Both': GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') else: - GroupWorkspaces(output_parameter_table_name, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') + GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_table_name, output_parameter_table_headers): @@ -614,8 +616,8 @@ def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_tabl return fit_workspace_name - def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers): - peak_params = self._find_peaks_and_output_params(spec_range, peak_position) + def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers): + peak_params = self._find_peaks_and_output_params(spec_range, peak_position, peak_index) fit_func = self._build_function_string(peak_params) start_str = 'composite=MultiDomainFunction, NumDeriv=1;' validSpecs = mtd[self._sample] @@ -630,7 +632,7 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam ties = ','.join(f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1,n_valid_specs)) func_string = composite_func + f';ties=({ties})' xmin, xmax = self._find_fit_x_window(peak_params) - fit_output_name = '__' + self._output_workspace_name + '_Peak_0' + fit_output_name = '__' + self._output_workspace_name + '_Peak_%d' % peak_index # create new workspace for each spectra x = validSpecs.readX(0) @@ -655,13 +657,17 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam self._output_fit_params_to_table_ws(0, fit_params, output_parameter_table_name, output_headers) + + fit_workspace_name = fws.name() self._del_fit_workspaces(ncm, fit_params, fws) - def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index=None, spec_index=None): - if spec_index is not None and peak_index is not None: + return fit_workspace_name + + def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index, spec_index=None): + if spec_index is not None: peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index) else: - peak_table_name = '__' + self._sample + '_peak_table' + peak_table_name = '__' + self._sample + '_peak_table_%d' % (peak_index) find_peak_params = self._get_find_peak_parameters(spec_number, [peak_position]) FindPeaks(InputWorkspace=self._sample, WorkspaceIndex=spec_index, PeaksList=peak_table_name, **find_peak_params) if mtd[peak_table_name].rowCount() == 0: diff --git a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py index 45e8aed0..b742ec5a 100644 --- a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py +++ b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py @@ -211,9 +211,27 @@ def test_copper_with_individual_and_global_fit(self, load_file_mock): detector_specific_r_tols = {"L1": {116: 0.65, 170: 0.75}} detector_specific_r_tols["L1"].update({k: TestConstants.INVALID_DETECTOR for k in [138, 141, 146, 147, 156, 158, 160, 163, 164, 165, 167, 168, 169, 170, 182, 191, 192]}) - self._assert_parameters_match_expected(params_table, detector_specific_r_tols) + self._assert_parameters_match_expected(params_table, detector_specific_r_tols, "Both") self._assert_E1_parameters_match_expected(params_table, 3e-3, "Both") + @patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._load_file') + def test_copper_with_global_fit(self, load_file_mock): + self._setup_copper_test() + self._output_workspace = "copper_analysis_test" + + load_file_mock.side_effect = self._load_file_side_effect + + self._create_evs_calibration_alg() + self._alg.setProperty("SharedParameterFitType", "Shared") + params_table = self._execute_evs_calibration_analysis() + + # Specify detectors tolerances set by user, then update with those to mask as invalid. + detector_specific_r_tols = {"L1": {116: 0.65, 170: 0.75}} + detector_specific_r_tols["L1"].update({k: TestConstants.INVALID_DETECTOR for k in [138, 141, 146, 147, 156, 158, 160, 163, 164, + 165, 167, 168, 169, 170, 182, 191, 192]}) + self._assert_parameters_match_expected(params_table, detector_specific_r_tols, "Shared") + self._assert_E1_parameters_match_expected(params_table, 3e-3, "Shared") + def _assert_theta_parameters_match_expected(self, params_table, rel_tolerance): thetas = params_table.column('theta') actual_thetas = self._calibrated_params['theta'] @@ -246,10 +264,12 @@ def _assert_E1_parameters_match_expected(self, params_table, rel_tolerance, fit_ global_E1 = params_table.column('global_E1')[0] self.assertAlmostEqual(global_E1, EVSGlobals.ENERGY_ESTIMATE, delta=EVSGlobals.ENERGY_ESTIMATE*rel_tolerance) - def _assert_parameters_match_expected(self, params_table, tolerances=None): + def _assert_parameters_match_expected(self, params_table, tolerances=None, fit_type="Individual"): rel_tol_theta, rel_tol_L1 = self._extract_tolerances(deepcopy(tolerances)) theta_errors = self._assert_theta_parameters_match_expected(params_table, rel_tol_theta) - L1_errors = self._assert_L1_parameters_match_expected(params_table, rel_tol_L1) + L1_errors = None + if fit_type != 'Shared': + L1_errors = self._assert_L1_parameters_match_expected(params_table, rel_tol_L1) if theta_errors or L1_errors: raise AssertionError(f"Theta: {theta_errors})\n L1: {L1_errors}") diff --git a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py index 04a7d572..6591a394 100644 --- a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py @@ -660,29 +660,32 @@ def test_fit_peaks_individual(self, group_workspaces_mock): group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters') - @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak') @patch('calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces') - def test_fit_peaks_shared(self, group_workspaces_mock, fit_shared_peak_mock): + def test_fit_peaks_shared(self, group_workspaces_mock): alg = EVSCalibrationFit() alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]])) alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) - alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') + alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}') alg._output_workspace_name = 'output_ws_name' alg._shared_parameter_fit_type = 'Shared' alg._spec_list = [3,4,5] alg._fit_peaks() - alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Shared_Peak_Parameters')]) - group_workspaces_mock.assert_called_once_with('output_ws_name_Shared_Peak_Parameters', + self.assertEqual(['fit_ws_0', 'fit_ws_1'], + alg._peak_fit_workspaces) + alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Shared_Peak_0_Parameters'), + call('output_ws_name_Shared_Peak_1_Parameters')]) + group_workspaces_mock.assert_called_once_with(['output_ws_name_Shared_Peak_0_Parameters','output_ws_name_Shared_Peak_1_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters') - fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c']) + alg._fit_shared_peak.assert_has_calls([call(0, 3, 10, 'output_ws_name_Shared_Peak_0_Parameters', ['a', 'b', 'c']), + call(1, 3, 41/6, 'output_ws_name_Shared_Peak_1_Parameters', ['a', 'b', 'c'])]) - @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak') @patch('calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces') - def test_fit_peaks_both(self, group_workspaces_mock, fit_shared_peak_mock): + def test_fit_peaks_both(self, group_workspaces_mock): alg = EVSCalibrationFit() alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]])) alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') + alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}') alg._output_workspace_name = 'output_ws_name' alg._shared_parameter_fit_type = 'Both' alg._spec_list = [3,4,5] @@ -695,11 +698,13 @@ def test_fit_peaks_both(self, group_workspaces_mock, fit_shared_peak_mock): call(1, 0, 2.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']), call(1, 1, 7.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']), call(1, 2, 10.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c'])]) - self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2']], + self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2'], 'fit_ws_0', 'fit_ws_1'], alg._peak_fit_workspaces) group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters', - 'output_ws_name_Shared_Peak_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters') - fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c']) + 'output_ws_name_Shared_Peak_0_Parameters','output_ws_name_Shared_Peak_1_Parameters'], + OutputWorkspace='output_ws_name_Peak_Parameters') + alg._fit_shared_peak.assert_has_calls([call(0, 3, 10, 'output_ws_name_Shared_Peak_0_Parameters', ['a', 'b', 'c']), + call(1, 3, 41/6, 'output_ws_name_Shared_Peak_1_Parameters', ['a', 'b', 'c'])]) @staticmethod def _setup_alg_mocks_fit_peak(): From 76510f508c6218d9904fbf998d275658433f24a2 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Thu, 4 May 2023 14:55:00 +0000 Subject: [PATCH 06/11] Exclude invalid spectra from global fit if preset --- .../calibration_scripts/calibrate_vesuvio_fit.py | 11 +++++++---- .../tests/system/test_system_analysis.py | 2 +- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index bc534967..d77d15a7 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -620,10 +620,13 @@ def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_paramet peak_params = self._find_peaks_and_output_params(spec_range, peak_position, peak_index) fit_func = self._build_function_string(peak_params) start_str = 'composite=MultiDomainFunction, NumDeriv=1;' - validSpecs = mtd[self._sample] - # need to decide how to remove invalid spectra from multifit - # MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) - # validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') + sample_ws = mtd[self._sample] + if self._invalid_detectors._detectors_preset == True: + invalid_spectra = self._invalid_detectors.get_invalid_detectors_index(self._spec_list) + MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) + validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') + else: + validSpecs = sample_ws n_valid_specs = validSpecs.getNumberHistograms() fit_func = ('(composite=CompositeFunction, NumDeriv=false, $domains=i;' + fit_func[:-1] + ');') * n_valid_specs diff --git a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py index b742ec5a..b119ec65 100644 --- a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py +++ b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py @@ -214,7 +214,7 @@ def test_copper_with_individual_and_global_fit(self, load_file_mock): self._assert_parameters_match_expected(params_table, detector_specific_r_tols, "Both") self._assert_E1_parameters_match_expected(params_table, 3e-3, "Both") - @patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._load_file') + @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._load_file') def test_copper_with_global_fit(self, load_file_mock): self._setup_copper_test() self._output_workspace = "copper_analysis_test" From 94f6aba91ce7b8dd5a8fb40e375137b13477f3c2 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Tue, 23 May 2023 14:36:22 +0000 Subject: [PATCH 07/11] address review comments --- .../calibrate_vesuvio_analysis.py | 174 +++++++++--------- .../calibrate_vesuvio_fit.py | 80 ++++---- 2 files changed, 125 insertions(+), 129 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py index f4259e56..381ad90a 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py @@ -250,56 +250,73 @@ def _calculate_final_flight_path(self, peak_table, spec_list, fit_type): Calculate the final flight path using the values for energy. This also uses the old value for L1 loaded from the parameter file. + @param peak_table - names of tables containing fitted parameters @param spec_list - spectrum range to calculate t0 for. + @param fit_type - type of fitting, one of 'Individual', 'Shared', or 'Both' """ if fit_type == 'Both': indv_peak_table = peak_table[0] - glob_peak_table = peak_table[1] - E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list) - E1 *= EVSGlobals.MEV_CONVERSION - global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list) - global_E1 *= EVSGlobals.MEV_CONVERSION + global_peak_table = peak_table[1] + params = self._read_calibration_params_from_table(spec_list, E1_col_names=['E1', 'global_E1']) elif fit_type == 'Shared': - glob_peak_table = peak_table[0] - global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list) - global_E1 *= EVSGlobals.MEV_CONVERSION + global_peak_table = peak_table[0] + params = self._read_calibration_params_from_table(spec_list, E1_col_names=['global_E1']) else: indv_peak_table = peak_table[0] - E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list) - E1 *= EVSGlobals.MEV_CONVERSION - - t0 = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list) - t0_error = EVSMiscFunctions.read_table_column(self._current_workspace, 't0_Err', spec_list) - L0 = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list) - theta = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list) - r_theta = EVSMiscFunctions.calculate_r_theta(self._sample_mass, theta) + params = self._read_calibration_params_from_table(spec_list, E1_col_names=['E1']) if fit_type != 'Shared': - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) - - delta_t = (peak_centres - t0) / 1e+6 - delta_t_error = t0_error / 1e+6 - v1 = np.sqrt( 2 * E1 /scipy.constants.m_n) - L1 = v1 * delta_t - L0 * r_theta - L1_error = v1 * delta_t_error - - self._set_table_column(self._current_workspace, 'L1', L1, spec_list) - self._set_table_column(self._current_workspace, 'L1_Err', L1_error, spec_list) + self._calculate_L1_from_table(spec_list, peak_centres, params, 'L1') if fit_type != 'Individual': - peak_centre = EVSMiscFunctions.read_fitting_result_table_column(glob_peak_table, 'f1.LorentzPos', [0]) - peak_centre_list = np.empty(spec_list[1] + 1 -spec_list[0]) + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(global_peak_table, 'f1.LorentzPos', [0]) + peak_centre_list = np.empty(spec_list[1] + 1 - spec_list[0]) peak_centre_list.fill(float(peak_centre)) + self._calculate_L1_from_table(spec_list, peak_centre_list, params, 'global_L1') + + def _read_calibration_params_from_table(self, spec_list, E1_col_names=[], L1_col_names=[]): + """ + Retrieve calibration parameters and return dict of param values + + @param spec_list - spectrum range to retreive param values for + @param E1_col_names - E1 columns to read from + @param L1_col_names - L1 columns to read from + """ + param_dict = {} + for col_name in E1_col_names: + if col_name == 'E1': + param_dict['E1'] = EVSMiscFunctions.read_table_column(self._current_workspace, col_name, spec_list) * EVSGlobals.MEV_CONVERSION + if col_name == 'global_E1': + param_dict['global_E1'] = EVSMiscFunctions.read_table_column(self._current_workspace, col_name, spec_list) * EVSGlobals.MEV_CONVERSION + for col_name in L1_col_names: + if col_name == 'L1': + param_dict['L1'] = EVSMiscFunctions.read_table_column(self._param_table, col_name, spec_list) + param_dict['t0'] = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list) + param_dict['t0_error'] = EVSMiscFunctions.read_table_column(self._current_workspace, 't0_Err', spec_list) + param_dict['L0'] = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list) + param_dict['theta'] = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list) + param_dict['r_theta'] = EVSMiscFunctions.calculate_r_theta(self._sample_mass, param_dict['theta']) + + return param_dict + + def _calculate_L1_from_table(self, spec_list, peak_centres, params, L1_col_name): + """ + Calculate L1 from fitted calibration parameters - delta_t = (peak_centre_list - t0) / 1e+6 - delta_t_error = t0_error / 1e+6 - global_v1 = np.sqrt( 2 * global_E1 /scipy.constants.m_n) - global_L1 = global_v1 * delta_t - L0 * r_theta - global_L1_error = global_v1 * delta_t_error + @param spec_list - spectrum range to calculate L1 for + @param peak_centres - list of peak centre values + @param params - dict of params required for calculation + @param L1_col_name - either 'L1' or 'global_L1' + """ + delta_t = (peak_centres - params['t0']) / 1e+6 + delta_t_error = params['t0_error'] / 1e+6 + v1 = np.sqrt( 2 * params[L1_col_name.strip('L1') + 'E1'] /scipy.constants.m_n) + L1 = v1 * delta_t - params['L0'] * params['r_theta'] + L1_error = v1 * delta_t_error - self._set_table_column(self._current_workspace, 'global_L1', global_L1, spec_list) - self._set_table_column(self._current_workspace, 'global_L1_Err', global_L1_error, spec_list) + self._set_table_column(self._current_workspace, L1_col_name, L1, spec_list) + self._set_table_column(self._current_workspace, f'{L1_col_name}_Err', L1_error, spec_list) def _calculate_scattering_angle(self, table_name, spec_list): """ @@ -347,70 +364,30 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): """ Calculate the final energy using the fitted peak centres of a run. - @param table_name - name of table containing fitted parameters for the peak centres + @param peak_table - name of table containing fitted parameters for the peak centres @param spec_list - spectrum range to calculate t0 for. + @param fit_type - type of fitting, one of 'Individual', 'Shared', or 'Both' """ - spec_range = EVSGlobals.DETECTOR_RANGE[1] + 1 - EVSGlobals.DETECTOR_RANGE[0] if not self._E1_value_and_error: - if fit_type == 'Both': indv_peak_table = peak_table[0] - glob_peak_table = peak_table[1] + global_peak_table = peak_table[1] elif fit_type == 'Shared': - glob_peak_table = peak_table[0] + global_peak_table = peak_table[0] else: indv_peak_table = peak_table[0] - - t0 = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list) - L0 = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list) - - L1 = EVSMiscFunctions.read_table_column(self._param_table, 'L1', spec_list) - theta = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list) - r_theta = EVSMiscFunctions.calculate_r_theta(self._sample_mass, theta) + + params = self._read_calibration_params_from_table(spec_list, L1_col_names=['L1']) if fit_type != 'Shared': - mean_E1 = np.empty(spec_range) - E1_error = np.empty(spec_range) - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) - - delta_t = (peak_centres - t0) / 1e+6 - v1 = (L0 * r_theta + L1) / delta_t - - E1 = 0.5*scipy.constants.m_n*v1**2 - E1 /= EVSGlobals.MEV_CONVERSION - - mean_E1_val = np.nanmean(E1) - E1_error_val = np.nanstd(E1) - - mean_E1.fill(mean_E1_val) - E1_error.fill(E1_error_val) - - self._set_table_column(self._current_workspace, 'E1', mean_E1) - self._set_table_column(self._current_workspace, 'E1_Err', E1_error) - + self._calculate_E1(spec_list, spec_range, peak_centres, params, 'E1') if fit_type != 'Individual': - global_E1 = np.empty(spec_range) - global_E1_error = np.empty(spec_range) - - peak_centre = EVSMiscFunctions.read_fitting_result_table_column(glob_peak_table, 'f1.LorentzPos', [0]) - peak_centre = [peak_centre] * spec_range - - delta_t = (peak_centre - t0) / 1e+6 - v1 = (L0 * r_theta + L1) / delta_t - E1 = 0.5*scipy.constants.m_n*v1**2 - E1 /= EVSGlobals.MEV_CONVERSION - - global_E1_val = np.nanmean(E1) - global_E1_error_val = np.nanstd(E1) - - global_E1.fill(global_E1_val) - global_E1_error.fill(global_E1_error_val) - - self._set_table_column(self._current_workspace, 'global_E1', global_E1) - self._set_table_column(self._current_workspace, 'global_E1_Err', global_E1_error) + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(global_peak_table, 'f1.LorentzPos', [0]) + peak_centre_list = [float(peak_centre)] * len(params['L0']) + self._calculate_E1(spec_list, spec_range, peak_centre_list, params, 'global_E1') else: mean_E1 = np.empty(spec_range) E1_error = np.empty(spec_range) @@ -424,6 +401,33 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): self._set_table_column(self._current_workspace, 'E1', mean_E1) self._set_table_column(self._current_workspace, 'E1_Err', E1_error) + def _calculate_E1(self, spec_list, spec_range, peak_centres, params, E1_col_name): + """ + Calculate E1 value from fitted parameters + + @param spec_list - spectrum range to calculate E1 for + @param spec_range - detector range + @param peak_centres - list of fitted peak centres + @param params - dict of params needed to calculate E1 + @param E1_col_name - either 'E1' or 'global_E1' + """ + E1_col = np.empty(spec_range) + E1_col_error = np.empty(spec_range) + + delta_t = (peak_centres - params['t0']) / 1e+6 + v1 = (params['L0'] * params['r_theta'] + params['L1']) / delta_t + E1 = 0.5 * scipy.constants.m_n * v1 ** 2 + E1 /= EVSGlobals.MEV_CONVERSION + + E1_val = np.nanmean(E1) + E1_error_val = np.nanstd(E1) + + E1_col.fill(E1_val) + E1_col_error.fill(E1_error_val) + + self._set_table_column(self._current_workspace, E1_col_name, E1_col) + self._set_table_column(self._current_workspace, f'{E1_col_name}_Err', E1_col_error) + def _setup(self): """ Setup algorithm. diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index d77d15a7..e129ca5f 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -574,25 +574,19 @@ def _fit_peaks(self): self._output_parameter_tables.append(output_parameter_table_name) self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) - if self._shared_parameter_fit_type == 'Individual': - GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') - if self._shared_parameter_fit_type != "Individual": for peak_index, estimated_peak_position in enumerate(np.mean(estimated_peak_positions_all_peaks,1)): - output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_%d_Parameters' % peak_index - output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name) + output_shared_parameter_table_name = self._output_workspace_name + '_Shared_Peak_%d_Parameters' % peak_index + output_shared_parameter_table_headers = self._create_parameter_table_and_output_headers(output_shared_parameter_table_name) - fit_workspace_name = self._fit_shared_peak(peak_index, self._spec_list[0], estimated_peak_position, output_parameter_table_name, - output_parameter_table_headers) + shared_fit_workspace_name = self._fit_shared_peak(peak_index, self._spec_list[0], estimated_peak_position, output_shared_parameter_table_name, + output_shared_parameter_table_headers) - self._output_parameter_tables.append(output_parameter_table_name) - self._peak_fit_workspaces.append(fit_workspace_name) + self._output_parameter_tables.append(output_shared_parameter_table_name) + self._peak_fit_workspaces.append(shared_fit_workspace_name) - if self._shared_parameter_fit_type == 'Both': - GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') - else: - GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') + GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters') def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_table_name, output_parameter_table_headers): @@ -618,42 +612,24 @@ def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_tabl def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers): peak_params = self._find_peaks_and_output_params(spec_range, peak_position, peak_index) - fit_func = self._build_function_string(peak_params) - start_str = 'composite=MultiDomainFunction, NumDeriv=1;' + sample_ws = mtd[self._sample] - if self._invalid_detectors._detectors_preset == True: - invalid_spectra = self._invalid_detectors.get_invalid_detectors_index(self._spec_list) - MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) - validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') - else: - validSpecs = sample_ws + invalid_spectra = self._invalid_detectors.identify_and_set_invalid_detectors_from_range(self._spec_list, output_parameter_table_name) + MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) + validSpecs = ExtractUnmaskedSpectra(InputWorkspace=sample_ws, OutputWorkspace='valid_spectra') n_valid_specs = validSpecs.getNumberHistograms() - - fit_func = ('(composite=CompositeFunction, NumDeriv=false, $domains=i;' + fit_func[:-1] + ');') * n_valid_specs - composite_func = start_str + fit_func[:-1] - - ties = ','.join(f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1,n_valid_specs)) - func_string = composite_func + f';ties=({ties})' + + func_string = self._generate_multifit_func_string(peak_params, n_valid_specs) + xmin, xmax = self._find_fit_x_window(peak_params) fit_output_name = '__' + self._output_workspace_name + '_Peak_%d' % peak_index - # create new workspace for each spectra - x = validSpecs.readX(0) - y = validSpecs.readY(0) - e = validSpecs.readE(0) - out_ws = self._sample + '_Spec_0' - CreateWorkspace(DataX=x, DataY=y, DataE=e, NSpec=1, OutputWorkspace=out_ws) - - other_inputs = [CreateWorkspace(DataX=validSpecs.readX(j), DataY=validSpecs.readY(j), DataE=validSpecs.readE(j), - OutputWorkspace=f'{self._sample}_Spec_{j}') - for j in range(1,n_valid_specs)] - - added_args = {f'InputWorkspace_{i + 1}': v for i,v in enumerate(other_inputs)} + multifit_workspaces = self._create_multifit_workspaces(validSpecs, n_valid_specs) - status, chi2, ncm, fit_params, fws, func, cost_func = Fit(Function=func_string, InputWorkspace=out_ws, IgnoreInvalidData=True, + status, chi2, ncm, fit_params, fws, func, cost_func = Fit(Function=func_string, InputWorkspace=f'{self._sample}_Spec_0', IgnoreInvalidData=True, StartX=xmin, EndX=xmax, CalcErrors=True, Output=fit_output_name, - Minimizer='Levenberg-Marquardt,RelError=1e-8', **added_args) + Minimizer='Levenberg-Marquardt,RelError=1e-8', **multifit_workspaces) [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0,n_valid_specs)] output_headers = ['f0.'+ name for name in output_parameter_table_headers] @@ -666,11 +642,27 @@ def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_paramet return fit_workspace_name + def _generate_multifit_func_string(self, peak_params, n_valid_specs): + start_str = 'composite=MultiDomainFunction, NumDeriv=1;' + fit_func = self._build_function_string(peak_params) + multifit_func = ('(composite=CompositeFunction, NumDeriv=false, $domains=i;' + fit_func[:-1] + ');') * n_valid_specs + ties = ','.join(f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1,n_valid_specs)) + func_string = start_str + multifit_func[:-1] + f';ties=({ties})' + return func_string + + def _create_multifit_workspaces(self, validSpecs, n_valid_specs): + multifit_inputs = [CreateWorkspace(DataX=validSpecs.readX(j), DataY=validSpecs.readY(j), DataE=validSpecs.readE(j), + OutputWorkspace=f'{self._sample}_Spec_{j}') + for j in range(0,n_valid_specs)] + + multifit_workspaces = {f'InputWorkspace_{i + 1}': v for i,v in enumerate(multifit_inputs[1:])} + + return multifit_workspaces + def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index, spec_index=None): + peak_table_name = f'__{self._sample}_peaks_table_{peak_index}' if spec_index is not None: - peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index) - else: - peak_table_name = '__' + self._sample + '_peak_table_%d' % (peak_index) + peak_table_name = peak_table_name + f'_{spec_index}' find_peak_params = self._get_find_peak_parameters(spec_number, [peak_position]) FindPeaks(InputWorkspace=self._sample, WorkspaceIndex=spec_index, PeaksList=peak_table_name, **find_peak_params) if mtd[peak_table_name].rowCount() == 0: From 8a366afd913430fa0c60e81b5f1e975762af6df1 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Tue, 23 May 2023 15:00:02 +0000 Subject: [PATCH 08/11] remove need to define indv/global peak tables --- .../calibrate_vesuvio_analysis.py | 60 +++++++++---------- 1 file changed, 27 insertions(+), 33 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py index 381ad90a..f8a41371 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py @@ -255,25 +255,15 @@ def _calculate_final_flight_path(self, peak_table, spec_list, fit_type): @param fit_type - type of fitting, one of 'Individual', 'Shared', or 'Both' """ if fit_type == 'Both': - indv_peak_table = peak_table[0] - global_peak_table = peak_table[1] params = self._read_calibration_params_from_table(spec_list, E1_col_names=['E1', 'global_E1']) + self._calculate_L1_from_table(spec_list, peak_table[0], params, 'L1') + self._calculate_L1_from_table(spec_list, peak_table[1], params, 'global_L1') elif fit_type == 'Shared': - global_peak_table = peak_table[0] params = self._read_calibration_params_from_table(spec_list, E1_col_names=['global_E1']) + self._calculate_L1_from_table(spec_list, peak_table[0], params, 'global_L1') else: - indv_peak_table = peak_table[0] params = self._read_calibration_params_from_table(spec_list, E1_col_names=['E1']) - - if fit_type != 'Shared': - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) - self._calculate_L1_from_table(spec_list, peak_centres, params, 'L1') - - if fit_type != 'Individual': - peak_centre = EVSMiscFunctions.read_fitting_result_table_column(global_peak_table, 'f1.LorentzPos', [0]) - peak_centre_list = np.empty(spec_list[1] + 1 - spec_list[0]) - peak_centre_list.fill(float(peak_centre)) - self._calculate_L1_from_table(spec_list, peak_centre_list, params, 'global_L1') + self._calculate_L1_from_table(spec_list, peak_table[0], params, 'L1') def _read_calibration_params_from_table(self, spec_list, E1_col_names=[], L1_col_names=[]): """ @@ -300,15 +290,22 @@ def _read_calibration_params_from_table(self, spec_list, E1_col_names=[], L1_col return param_dict - def _calculate_L1_from_table(self, spec_list, peak_centres, params, L1_col_name): + def _calculate_L1_from_table(self, spec_list, peak_table, params, L1_col_name): """ Calculate L1 from fitted calibration parameters @param spec_list - spectrum range to calculate L1 for - @param peak_centres - list of peak centre values + @param peak_table - name of table to calculate L1 for @param params - dict of params required for calculation @param L1_col_name - either 'L1' or 'global_L1' """ + if L1_col_name == 'L1': + peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, peak_table) + elif L1_col_name == 'global_L1': + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(peak_table, 'f1.LorentzPos', [0]) + peak_centres = np.empty(spec_list[1] + 1 - spec_list[0]) + peak_centres.fill(float(peak_centre)) + delta_t = (peak_centres - params['t0']) / 1e+6 delta_t_error = params['t0_error'] / 1e+6 v1 = np.sqrt( 2 * params[L1_col_name.strip('L1') + 'E1'] /scipy.constants.m_n) @@ -364,30 +361,21 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): """ Calculate the final energy using the fitted peak centres of a run. - @param peak_table - name of table containing fitted parameters for the peak centres + @param peak_table - list of tables containing fitted parameters for the peak centres @param spec_list - spectrum range to calculate t0 for. @param fit_type - type of fitting, one of 'Individual', 'Shared', or 'Both' """ spec_range = EVSGlobals.DETECTOR_RANGE[1] + 1 - EVSGlobals.DETECTOR_RANGE[0] if not self._E1_value_and_error: + params = self._read_calibration_params_from_table(spec_list, L1_col_names=['L1']) if fit_type == 'Both': - indv_peak_table = peak_table[0] - global_peak_table = peak_table[1] + self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'E1') + self._calculate_E1(spec_list, spec_range, peak_table[1], params, 'global_E1') elif fit_type == 'Shared': - global_peak_table = peak_table[0] + self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'global_E1') else: - indv_peak_table = peak_table[0] - - params = self._read_calibration_params_from_table(spec_list, L1_col_names=['L1']) - - if fit_type != 'Shared': - peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, indv_peak_table) - self._calculate_E1(spec_list, spec_range, peak_centres, params, 'E1') - if fit_type != 'Individual': - peak_centre = EVSMiscFunctions.read_fitting_result_table_column(global_peak_table, 'f1.LorentzPos', [0]) - peak_centre_list = [float(peak_centre)] * len(params['L0']) - self._calculate_E1(spec_list, spec_range, peak_centre_list, params, 'global_E1') + self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'E1') else: mean_E1 = np.empty(spec_range) E1_error = np.empty(spec_range) @@ -401,16 +389,22 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): self._set_table_column(self._current_workspace, 'E1', mean_E1) self._set_table_column(self._current_workspace, 'E1_Err', E1_error) - def _calculate_E1(self, spec_list, spec_range, peak_centres, params, E1_col_name): + def _calculate_E1(self, spec_list, spec_range, peak_table, params, E1_col_name): """ Calculate E1 value from fitted parameters @param spec_list - spectrum range to calculate E1 for @param spec_range - detector range - @param peak_centres - list of fitted peak centres + @param peak_table - name of table containing fitted parameters for the peak centres @param params - dict of params needed to calculate E1 @param E1_col_name - either 'E1' or 'global_E1' """ + if E1_col_name == 'E1': + peak_centres = self._invalid_detectors.filter_peak_centres_for_invalid_detectors(spec_list, peak_table) + elif E1_col_name == 'global_E1': + peak_centre = EVSMiscFunctions.read_fitting_result_table_column(peak_table, 'f1.LorentzPos', [0]) + peak_centres = [float(peak_centre)] * len(params['L0']) + E1_col = np.empty(spec_range) E1_col_error = np.empty(spec_range) From 41b42e1b74086b8b72e584263f296d3b7a58cbb0 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 24 May 2023 08:58:38 +0000 Subject: [PATCH 09/11] Compare enums instead of strings --- .../calibrate_vesuvio_analysis.py | 18 +++++++++--------- .../calibrate_vesuvio_fit.py | 10 +++++----- .../calibrate_vesuvio_helper_functions.py | 6 ++++++ 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py index f8a41371..d46fad35 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py @@ -4,7 +4,7 @@ from mantid.simpleapi import CreateEmptyTableWorkspace, DeleteWorkspace, ReplaceSpecialValues, GroupWorkspaces, mtd,\ ConvertTableToMatrixWorkspace, ConjoinWorkspaces, Transpose, PlotPeakByLogValue,RenameWorkspace from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, EVSMiscFunctions,\ - InvalidDetectors + InvalidDetectors, FitTypes import os import sys @@ -53,7 +53,7 @@ def PyInit(self): self.declareProperty('Iterations', 2, validator=IntBoundedValidator(lower=1), doc="Number of iterations to perform. Default is 2.") - shared_fit_type_validator = StringListValidator(["Individual", "Shared", "Both"]) + shared_fit_type_validator = StringListValidator([member.value for member in FitTypes]) self.declareProperty('SharedParameterFitType', "Individual", doc='Calculate shared parameters using an individual and/or' 'global fit.', validator=shared_fit_type_validator) @@ -132,7 +132,7 @@ def PyExec(self): self._run_calibration_fit(Samples=self._samples, Function='Voigt', Mode='SingleDifference', SpectrumRange=EVSGlobals.BACKSCATTERING_RANGE, InstrumentParameterWorkspace=self._param_table, Mass=self._sample_mass, OutputWorkspace=E1_fit_back, CreateOutput=self._create_output, - PeakType='Recoil', SharedParameterFitType=self._shared_parameter_fit_type, + PeakType='Recoil', SharedParameterFitType=self._shared_parameter_fit_type.value, InvalidDetectors=self._invalid_detectors.get_all_invalid_detectors()) self._invalid_detectors.add_invalid_detectors(invalid_detectors) @@ -148,7 +148,7 @@ def PyExec(self): self._run_calibration_fit(Samples=self._samples, Function='Voigt', Mode='SingleDifference', SpectrumRange=EVSGlobals.FRONTSCATTERING_RANGE, InstrumentParameterWorkspace=self._param_table, Mass=self._sample_mass, OutputWorkspace=E1_fit_front, CreateOutput=self._create_output, - PeakType='Recoil', SharedParameterFitType=self._shared_parameter_fit_type, + PeakType='Recoil', SharedParameterFitType=self._shared_parameter_fit_type.value, InvalidDetectors=self._invalid_detectors.get_all_invalid_detectors()) self._invalid_detectors.add_invalid_detectors(invalid_detectors) @@ -254,11 +254,11 @@ def _calculate_final_flight_path(self, peak_table, spec_list, fit_type): @param spec_list - spectrum range to calculate t0 for. @param fit_type - type of fitting, one of 'Individual', 'Shared', or 'Both' """ - if fit_type == 'Both': + if fit_type is FitTypes.BOTH: params = self._read_calibration_params_from_table(spec_list, E1_col_names=['E1', 'global_E1']) self._calculate_L1_from_table(spec_list, peak_table[0], params, 'L1') self._calculate_L1_from_table(spec_list, peak_table[1], params, 'global_L1') - elif fit_type == 'Shared': + elif fit_type is FitTypes.SHARED: params = self._read_calibration_params_from_table(spec_list, E1_col_names=['global_E1']) self._calculate_L1_from_table(spec_list, peak_table[0], params, 'global_L1') else: @@ -369,10 +369,10 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type): if not self._E1_value_and_error: params = self._read_calibration_params_from_table(spec_list, L1_col_names=['L1']) - if fit_type == 'Both': + if fit_type is FitTypes.BOTH: self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'E1') self._calculate_E1(spec_list, spec_range, peak_table[1], params, 'global_E1') - elif fit_type == 'Shared': + elif fit_type is FitTypes.SHARED: self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'global_E1') else: self._calculate_E1(spec_list, spec_range, peak_table[0], params, 'E1') @@ -433,7 +433,7 @@ def _setup(self): self._d_spacings = self.getProperty("DSpacings").value.tolist() self._E1_value_and_error = self.getProperty("E1FixedValueAndError").value.tolist() self._invalid_detectors = InvalidDetectors(self.getProperty("InvalidDetectors").value.tolist()) - self._shared_parameter_fit_type = self.getProperty("SharedParameterFitType").value + self._shared_parameter_fit_type = FitTypes(self.getProperty("SharedParameterFitType").value) self._calc_L0 = self.getProperty("CalculateL0").value self._make_IP_file = self.getProperty("CreateIPFile").value self._output_workspace_name = self.getPropertyValue("OutputWorkspace") diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index e129ca5f..3b2663ba 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -6,7 +6,7 @@ ReplaceSpecialValues, FindPeaks, GroupWorkspaces, mtd, Plus, LoadVesuvio, LoadRaw, ConvertToDistribution, FindPeakBackground,\ ExtractSingleSpectrum, SumSpectra, AppendSpectra, CloneWorkspace, Fit, MaskDetectors, ExtractUnmaskedSpectra, CreateWorkspace from functools import partial -from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, EVSMiscFunctions, InvalidDetectors +from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, EVSMiscFunctions, InvalidDetectors, FitTypes import os import sys @@ -63,7 +63,7 @@ def PyInit(self): doc='Choose the peak type that is being fitted.' 'Note that supplying a set of dspacings overrides the setting here') - shared_fit_type_validator = StringListValidator(["Individual", "Shared", "Both"]) + shared_fit_type_validator = StringListValidator([member.value for member in FitTypes]) self.declareProperty('SharedParameterFitType', "Individual", doc='Calculate shared parameters using an individual and/or' 'global fit.', validator=shared_fit_type_validator) @@ -122,7 +122,7 @@ def _setup_class_variables_from_properties(self): self._energy_estimates = self.getProperty("Energy").value self._sample_mass = self.getProperty("Mass").value self._create_output = self.getProperty("CreateOutput").value - self._shared_parameter_fit_type = self.getProperty("SharedParameterFitType").value + self._shared_parameter_fit_type = FitTypes(self.getProperty("SharedParameterFitType").value) self._invalid_detectors = InvalidDetectors(self.getProperty("InvalidDetectors").value.tolist()) def _setup_spectra_list(self): @@ -558,7 +558,7 @@ def _fit_peaks(self): num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape self._output_parameter_tables = [] self._peak_fit_workspaces = [] - if self._shared_parameter_fit_type != "Shared": + if self._shared_parameter_fit_type is not FitTypes.SHARED: self._prog_reporter = Progress(self, 0.0, 1.0, num_estimated_peaks*num_spectra) for peak_index, estimated_peak_positions in enumerate(estimated_peak_positions_all_peaks): @@ -574,7 +574,7 @@ def _fit_peaks(self): self._output_parameter_tables.append(output_parameter_table_name) self._peak_fit_workspaces.append(self._peak_fit_workspaces_by_spec) - if self._shared_parameter_fit_type != "Individual": + if self._shared_parameter_fit_type is not FitTypes.INDIVIDUAL: for peak_index, estimated_peak_position in enumerate(np.mean(estimated_peak_positions_all_peaks,1)): output_shared_parameter_table_name = self._output_workspace_name + '_Shared_Peak_%d_Parameters' % peak_index diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_helper_functions.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_helper_functions.py index ae86398e..af82c278 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_helper_functions.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_helper_functions.py @@ -1,5 +1,6 @@ import numpy as np import scipy.constants +from enum import Enum from mantid.simpleapi import CreateEmptyTableWorkspace, mtd @@ -73,6 +74,11 @@ class EVSGlobals: # 1 meV in Joules MEV_CONVERSION = 1.602176487e-22 +class FitTypes(Enum): + + INDIVIDUAL = "Individual" + SHARED = "Shared" + BOTH = "Both" class EVSMiscFunctions: From 0fb2ab8f9b9d1cc4ddde46e30cd8435c722a263f Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 24 May 2023 09:09:50 +0000 Subject: [PATCH 10/11] fix failing tests --- .../tests/unit/test_calibrate_vesuvio_fit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py index 6591a394..56a03522 100644 --- a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py @@ -1,5 +1,5 @@ from calibration_scripts.calibrate_vesuvio_fit import EVSCalibrationFit -from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals +from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, FitTypes from mock import MagicMock, patch, call from functools import partial from mantid.kernel import IntArrayProperty, StringArrayProperty, FloatArrayProperty @@ -645,7 +645,7 @@ def test_fit_peaks_individual(self, group_workspaces_mock): alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') alg._output_workspace_name = 'output_ws_name' - alg._shared_parameter_fit_type = 'Individual' + alg._shared_parameter_fit_type = FitTypes.INDIVIDUAL alg._fit_peaks() alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Peak_0_Parameters'), call('output_ws_name_Peak_1_Parameters')]) @@ -667,7 +667,7 @@ def test_fit_peaks_shared(self, group_workspaces_mock): alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c']) alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}') alg._output_workspace_name = 'output_ws_name' - alg._shared_parameter_fit_type = 'Shared' + alg._shared_parameter_fit_type = FitTypes.SHARED alg._spec_list = [3,4,5] alg._fit_peaks() self.assertEqual(['fit_ws_0', 'fit_ws_1'], @@ -687,7 +687,7 @@ def test_fit_peaks_both(self, group_workspaces_mock): alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}') alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}') alg._output_workspace_name = 'output_ws_name' - alg._shared_parameter_fit_type = 'Both' + alg._shared_parameter_fit_type = FitTypes.BOTH alg._spec_list = [3,4,5] alg._fit_peaks() alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Peak_0_Parameters'), From 5f4ef97b3220dc81ff272824f3c9ebef010adf79 Mon Sep 17 00:00:00 2001 From: Jess Farmer Date: Wed, 31 May 2023 08:32:19 +0000 Subject: [PATCH 11/11] add more unit tests for _fit_shared_peak --- .../calibrate_vesuvio_fit.py | 10 +-- .../tests/unit/test_calibrate_vesuvio_fit.py | 87 ++++++++++++++++++- 2 files changed, 91 insertions(+), 6 deletions(-) diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py index 3b2663ba..16024a0b 100644 --- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py @@ -612,7 +612,7 @@ def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_tabl def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers): peak_params = self._find_peaks_and_output_params(spec_range, peak_position, peak_index) - + sample_ws = mtd[self._sample] invalid_spectra = self._invalid_detectors.identify_and_set_invalid_detectors_from_range(self._spec_list, output_parameter_table_name) MaskDetectors(Workspace=sample_ws, WorkspaceIndexList=invalid_spectra) @@ -630,15 +630,13 @@ def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_paramet StartX=xmin, EndX=xmax, CalcErrors=True, Output=fit_output_name, Minimizer='Levenberg-Marquardt,RelError=1e-8', **multifit_workspaces) - [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0,n_valid_specs)] - output_headers = ['f0.'+ name for name in output_parameter_table_headers] self._output_fit_params_to_table_ws(0, fit_params, output_parameter_table_name, output_headers) fit_workspace_name = fws.name() - self._del_fit_workspaces(ncm, fit_params, fws) + self._del_fit_workspaces(ncm, fit_params, fws, n_valid_specs) return fit_workspace_name @@ -693,11 +691,13 @@ def _output_fit_params_to_table_ws(self, spec_num, params, output_table_name, ou mtd[output_table_name].addRow([spec_num] + row) - def _del_fit_workspaces(self, ncm, fit_params, fws): + def _del_fit_workspaces(self, ncm, fit_params, fws, n_valid_specs=None): DeleteWorkspace(ncm) DeleteWorkspace(fit_params) if not self._create_output: DeleteWorkspace(fws) + if n_valid_specs != None: + [DeleteWorkspace(f"{self._sample}_Spec_{i}") for i in range(0,n_valid_specs)] def _get_find_peak_parameters(self, spec_number, peak_centre, unconstrained=False): """ diff --git a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py index 56a03522..ec9a3378 100644 --- a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py +++ b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py @@ -1,5 +1,5 @@ from calibration_scripts.calibrate_vesuvio_fit import EVSCalibrationFit -from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, FitTypes +from calibration_scripts.calibrate_vesuvio_helper_functions import EVSGlobals, FitTypes, InvalidDetectors from mock import MagicMock, patch, call from functools import partial from mantid.kernel import IntArrayProperty, StringArrayProperty, FloatArrayProperty @@ -745,6 +745,83 @@ def test_fit_peak(self, mock_fit): self.assertEqual(fit_workspace_name, 'fws') + @staticmethod + def _setup_alg_mocks_fit_shared_peak(): + alg = EVSCalibrationFit() + alg._find_peaks_and_output_params = MagicMock(return_value='peak_params') + alg._generate_multifit_func_string = MagicMock(return_value='function_string') + alg._find_fit_x_window = MagicMock(return_value=(0, 10)) + alg._create_multifit_workspaces = MagicMock(return_value={'InputWorkspace_1': 'workspace'}) + alg._output_fit_params_to_table_ws = MagicMock() + alg._sample = 'sample' + alg._output_workspace_name = 'output' + alg._spec_list = [30] + alg._prog_reporter = MagicMock() + alg._invalid_detectors = InvalidDetectors([]) + alg._del_fit_workspaces = MagicMock() + + return alg + + @patch('calibration_scripts.calibrate_vesuvio_fit.Fit') + @patch('calibration_scripts.calibrate_vesuvio_fit.InvalidDetectors.identify_and_set_invalid_detectors_from_range') + @patch('calibration_scripts.calibrate_vesuvio_fit.mtd') + @patch('calibration_scripts.calibrate_vesuvio_fit.MaskDetectors') + @patch('calibration_scripts.calibrate_vesuvio_fit.ExtractUnmaskedSpectra') + def test_fit_shared_peak(self, mock_extract, mock_mask, mock_mtd, mock_invalid, mock_fit): + alg = self._setup_alg_mocks_fit_shared_peak() + + mock_mtd.__getitem__.return_value = 'sample' + mock_invalid.return_value = [] + validSpecs = MagicMock() + validSpecs.getNumberHistograms.return_value = 1 + mock_extract.return_value = validSpecs + fit_workspace = MagicMock() + fit_workspace.name.return_value = 'fws' + fit_results = ('success', 0.5, 'ncm', 'fit_params', fit_workspace, 'func', 'cost_func') + mock_fit.return_value = fit_results + + fit_workspace_name = alg._fit_shared_peak(0, 30, 5, 'output_Peak_0_Parameters', ['a', 'b', 'c']) + + alg._find_peaks_and_output_params.assert_called_once_with(30, 5, 0) + mock_invalid.assert_called_once_with([30], 'output_Peak_0_Parameters') + mock_mask.assert_called_once_with(Workspace='sample', WorkspaceIndexList=[]) + mock_extract.assert_called_once_with(InputWorkspace='sample', OutputWorkspace='valid_spectra') + alg._generate_multifit_func_string.assert_called_once_with('peak_params', 1) + alg._find_fit_x_window.assert_called_once_with('peak_params') + mock_fit.assert_called_once_with(Function='function_string', InputWorkspace='sample_Spec_0', + IgnoreInvalidData=True, StartX=0, EndX=10, + CalcErrors=True, Output='__output_Peak_0', + Minimizer='Levenberg-Marquardt,RelError=1e-8', + InputWorkspace_1='workspace') + alg._output_fit_params_to_table_ws.assert_called_once_with(0, 'fit_params', 'output_Peak_0_Parameters', + ['f0.a', 'f0.b', 'f0.c']) + alg._del_fit_workspaces.assert_called_once_with('ncm', 'fit_params', fit_workspace, 1) + + self.assertEqual(fit_workspace_name, 'fws') + + @patch('calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._build_function_string') + def test_generate_multifit_func_string(self, mock_build_func_string): + alg = EVSCalibrationFit() + alg._func_param_names = {'key1': 'val1', 'key2': 'val2', 'key3': 'val3'} + mock_build_func_string.return_value = 'function_string;' + func_string = alg._generate_multifit_func_string('peak_params', 2) + expected = ("composite=MultiDomainFunction, NumDeriv=1;(composite=CompositeFunction, " + "NumDeriv=false, $domains=i;function_string);(composite=CompositeFunction, NumDeriv=false, " + "$domains=i;function_string);ties=(f1.f1.val1=f0.f1.val1,f1.f1.val2=f0.f1.val2,f1.f1.val3=f0.f1.val3)") + self.assertEqual(func_string, expected) + + + @patch('calibration_scripts.calibrate_vesuvio_fit.CreateWorkspace') + def test_create_multifit_workspaces(self, create_workspace_mock): + alg = EVSCalibrationFit() + alg._sample = 'sample' + validSpecs = MagicMock() + validSpecs.readX.return_value = 0 + validSpecs.readY.return_value = 0 + validSpecs.readE.return_value = 0 + alg._create_multifit_workspaces(validSpecs, 1) + create_workspace_mock.assert_has_calls([call(DataX=0,DataY=0,DataE=0, OutputWorkspace='sample_Spec_0')]) + @patch('calibration_scripts.calibrate_vesuvio_fit.mtd') @patch('calibration_scripts.calibrate_vesuvio_fit.DeleteWorkspace') @patch('calibration_scripts.calibrate_vesuvio_fit.FindPeaks') @@ -850,6 +927,14 @@ def test_del_fit_workspace_create_output_true(self, del_ws_mock): alg._del_fit_workspaces('ncm', 'params', 'fws') del_ws_mock.assert_has_calls([call('ncm'), call('params'), call('fws')]) + @patch('calibration_scripts.calibrate_vesuvio_fit.DeleteWorkspace') + def test_del_fit_workspace_shared_fit(self, del_ws_mock): + alg = EVSCalibrationFit() + alg._sample = 'sample' + alg._create_output = False + alg._del_fit_workspaces('ncm', 'params', 'fws', 1) + del_ws_mock.assert_has_calls([call('ncm'), call('params'), call('fws'), call('sample_Spec_0')]) + @patch('calibration_scripts.calibrate_vesuvio_fit.CreateEmptyTableWorkspace') @patch('calibration_scripts.calibrate_vesuvio_fit.AnalysisDataService') def test_create_output_parameters_table_ws(self, mock_ADS, mock_create_empty_table_ws):