diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py index 1dcbe0b..a68a77e 100755 --- a/CHAP/edd/processor.py +++ b/CHAP/edd/processor.py @@ -796,6 +796,7 @@ def get_nxroot(self, from CHAP.edd.utils import ( get_peak_locations, get_unique_hkls_ds, + get_spectra_fits ) from CHAP.utils.fit import FitMap @@ -976,52 +977,20 @@ def linkdims(nxgroup, field_dims=[]): fit_ds = np.asarray([ds[i] for i in detector.hkl_indices]) peak_locations = get_peak_locations( fit_ds, detector.tth_calibrated) - num_peak = len(peak_locations) - # KLS: Use the below def of peak_locations when - # FitMap.create_multipeak_model can accept a list of maps - # for centers. - # tth = np.radians(detector.map_tth(map_config)) - # peak_locations = [get_peak_locations(d0, tth) for d0 in fit_ds] - - # Perform initial fit: assume uniform strain for all HKLs - self.logger.debug('Performing uniform fit') - fit = FitMap(det_nxdata.intensity.nxdata, x=energies) - delta = 0.1 * (energies[-1]-energies[0]) - centers_range = ( - max(0.0, energies[0]-delta), energies[-1]+delta) - fit.create_multipeak_model( - peak_locations, - fit_type='uniform', - peak_models=detector.peak_models, - background=detector.background, - fwhm_min=detector.fwhm_min, - fwhm_max=detector.fwhm_max, - centers_range=centers_range) - fit.fit() - uniform_fit_centers = [ - fit.best_values[ - fit.best_parameters().index(f'peak{i+1}_center')] - for i in range(num_peak)] - uniform_fit_centers_errors = [ - fit.best_errors[ - fit.best_parameters().index(f'peak{i+1}_center')] - for i in range(num_peak)] - uniform_fit_amplitudes = [ - fit.best_values[ - fit.best_parameters().index(f'peak{i+1}_amplitude')] - for i in range(num_peak)] - uniform_fit_amplitudes_errors = [ - fit.best_errors[ - fit.best_parameters().index(f'peak{i+1}_amplitude')] - for i in range(num_peak)] - uniform_fit_sigmas = [ - fit.best_values[ - fit.best_parameters().index(f'peak{i+1}_sigma')] - for i in range(num_peak)] - uniform_fit_sigmas_errors = [ - fit.best_errors[ - fit.best_parameters().index(f'peak{i+1}_sigma')] - for i in range(num_peak)] + + (uniform_fit_centers, uniform_fit_centers_errors, + uniform_fit_amplitudes, uniform_fit_amplitudes_errors, + uniform_fit_sigmas, uniform_fit_sigmas_errors, + uniform_best_fit, uniform_residuals, + uniform_redchi, uniform_success, + unconstrained_fit_centers, unconstrained_fit_centers_errors, + unconstrained_fit_amplitudes, unconstrained_fit_amplitudes_errors, + unconstrained_fit_sigmas, unconstrained_fit_sigmas_errors, + unconstrained_best_fit, unconstrained_residuals, + unconstrained_redchi, unconstrained_success) = \ + get_spectra_fits( + det_nxdata.intensity.nxdata, energies, + peak_locations, detector) # Add uniform fit results to the NeXus structure nxdetector.uniform_fit = NXcollection() @@ -1033,10 +1002,10 @@ def linkdims(nxgroup, field_dims=[]): linkdims( fit_nxdata, {'axes': 'energy', 'index': len(map_config.shape)}) fit_nxdata.makelink(det_nxdata.energy) - fit_nxdata.best_fit= fit.best_fit - fit_nxdata.residuals = fit.residual - fit_nxdata.redchi = fit.redchi - fit_nxdata.success = fit.success + fit_nxdata.best_fit= uniform_best_fit + fit_nxdata.residuals = uniform_residuals + fit_nxdata.redchi = uniform_redchi + fit_nxdata.success = uniform_success # Peak-by-peak results # fit_nxgroup.fit_hkl_centers = NXdata() @@ -1082,39 +1051,6 @@ def linkdims(nxgroup, field_dims=[]): value=sigmas_error) fit_nxgroup[hkl_name].sigmas.attrs['signal'] = 'values' - # Perform second fit: do not assume uniform strain for all - # HKLs, and use the fit peak centers from the uniform fit - # as inital guesses - self.logger.debug('Performing unconstrained fit') - fit.create_multipeak_model(fit_type='unconstrained') - fit.fit(rel_amplitude_cutoff=detector.rel_amplitude_cutoff) - unconstrained_fit_centers = np.array( - [fit.best_values[ - fit.best_parameters()\ - .index(f'peak{i+1}_center')] - for i in range(num_peak)]) - unconstrained_fit_centers_errors = np.array( - [fit.best_errors[ - fit.best_parameters()\ - .index(f'peak{i+1}_center')] - for i in range(num_peak)]) - unconstrained_fit_amplitudes = [ - fit.best_values[ - fit.best_parameters().index(f'peak{i+1}_amplitude')] - for i in range(num_peak)] - unconstrained_fit_amplitudes_errors = [ - fit.best_errors[ - fit.best_parameters().index(f'peak{i+1}_amplitude')] - for i in range(num_peak)] - unconstrained_fit_sigmas = [ - fit.best_values[ - fit.best_parameters().index(f'peak{i+1}_sigma')] - for i in range(num_peak)] - unconstrained_fit_sigmas_errors = [ - fit.best_errors[ - fit.best_parameters().index(f'peak{i+1}_sigma')] - for i in range(num_peak)] - if interactive or save_figures: # Third party modules import matplotlib.animation as animation @@ -1131,9 +1067,10 @@ def animate(i): intensity.set_ydata( det_nxdata.intensity.nxdata[map_index] / det_nxdata.intensity.nxdata[map_index].max()) - best_fit.set_ydata(fit.best_fit[map_index] - / fit.best_fit[map_index].max()) - # residual.set_ydata(fit.residual[map_index]) + best_fit.set_ydata( + unconstrained_best_fit[map_index] + / unconstrained_best_fit[map_index].max()) + # residual.set_ydata(unconstrained_residuals[map_index]) index.set_text('\n'.join(f'{k}[{i}] = {v}' for k, v in map_config.get_coords(map_index).items())) if save_figures: @@ -1149,12 +1086,12 @@ def animate(i): / det_nxdata.intensity.nxdata[map_index].max()) intensity, = ax.plot( energies, data_normalized, 'b.', label='data') - fit_normalized = (fit.best_fit[map_index] - / fit.best_fit[map_index].max()) + fit_normalized = (unconstrained_best_fit[map_index] + / unconstrained_best_fit[map_index].max()) best_fit, = ax.plot( energies, fit_normalized, 'k-', label='fit') # residual, = ax.plot( - # energies, fit.residual[map_index], 'r-', + # energies, unconstrained_residuals[map_index], 'r-', # label='residual') ax.set( title='Unconstrained fits', @@ -1225,10 +1162,10 @@ def animate(i): linkdims( fit_nxdata, {'axes': 'energy', 'index': len(map_config.shape)}) fit_nxdata.makelink(det_nxdata.energy) - fit_nxdata.best_fit= fit.best_fit - fit_nxdata.residuals = fit.residual - fit_nxdata.redchi = fit.redchi - fit_nxdata.success = fit.success + fit_nxdata.best_fit= unconstrained_best_fit + fit_nxdata.residuals = unconstrained_residuals + fit_nxdata.redchi = unconstrained_redchi + fit_nxdata.success = unconstrained_success # Peak-by-peak results fit_nxgroup.fit_hkl_centers = NXdata() diff --git a/CHAP/edd/utils.py b/CHAP/edd/utils.py index b5d2c30..44345b6 100755 --- a/CHAP/edd/utils.py +++ b/CHAP/edd/utils.py @@ -873,3 +873,127 @@ def confirm(event): fig.tight_layout(rect=(0, 0, 1, 0.95)) return fig, selected_bin_ranges, selected_hkl_indices + + +def get_spectra_fits(spectra, energies, peak_locations, fit_params): + """Return twenty arrays of fit results for the map of spectra + provided: uniform centers, uniform center errors, uniform + amplitudes, uniform amplitude errors, uniform sigmas, uniform + sigma errors, uniform best fit, uniform residuals, uniform reduced + chi, uniform success codes, unconstrained centers, unconstrained + center errors, unconstrained amplitudes, unconstrained amplitude + errors, unconstrained sigmas, unconstrained sigma errors, + unconstrained best fit, unconstrained residuals, unconstrained + reduced chi, and unconstrained success codes. + + :param spectra: Array of intensity spectra to fit. + :type spectra: numpy.ndarray + :param energies: Bin energies for the spectra provided. + :type energies: numpy.ndarray + :param peak_locations: Initial guesses for peak ceneters to use + for the uniform fit. + :type peak_locations: list[float] + :param fit_params: Detector element fit parameters. + :type fit_params: CHAP.edd.models.MCAElementStrainAnalysisConfig + :returns: Uniform and unconstrained centers, amplitdues, sigmas + (and errors for all three), best fits, residuals between the + best fits and the input spectra, reduced chi, and fit success + statuses. + :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, + numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, + numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, + numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, + numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, + numpy.ndarray] + """ + from CHAP.utils.fit import FitMap + + # Perform fit to get measured peak positions + fit = FitMap(spectra, x=energies) + num_peak = len(peak_locations) + delta = 0.1 * (energies[-1]-energies[0]) + centers_range = ( + max(0.0, energies[0]-delta), energies[-1]+delta) + fit.create_multipeak_model( + peak_locations, + fit_type='uniform', + peak_models=fit_params.peak_models, + background=fit_params.background, + fwhm_min=fit_params.fwhm_min, + fwhm_max=fit_params.fwhm_max, + centers_range=centers_range) + fit.fit() + uniform_fit_centers = [ + fit.best_values[ + fit.best_parameters().index(f'peak{i+1}_center')] + for i in range(num_peak)] + uniform_fit_centers_errors = [ + fit.best_errors[ + fit.best_parameters().index(f'peak{i+1}_center')] + for i in range(num_peak)] + uniform_fit_amplitudes = [ + fit.best_values[ + fit.best_parameters().index(f'peak{i+1}_amplitude')] + for i in range(num_peak)] + uniform_fit_amplitudes_errors = [ + fit.best_errors[ + fit.best_parameters().index(f'peak{i+1}_amplitude')] + for i in range(num_peak)] + uniform_fit_sigmas = [ + fit.best_values[ + fit.best_parameters().index(f'peak{i+1}_sigma')] + for i in range(num_peak)] + uniform_fit_sigmas_errors = [ + fit.best_errors[ + fit.best_parameters().index(f'peak{i+1}_sigma')] + for i in range(num_peak)] + uniform_best_fit = fit.best_fit + uniform_residuals = fit.residual + uniform_redchi = fit.redchi + uniform_success = fit.success + + # Perform unconstrained fit + fit.create_multipeak_model(fit_type='unconstrained') + fit.fit(rel_amplitude_cutoff=fit_params.rel_amplitude_cutoff) + unconstrained_fit_centers = np.array( + [fit.best_values[ + fit.best_parameters()\ + .index(f'peak{i+1}_center')] + for i in range(num_peak)]) + unconstrained_fit_centers_errors = np.array( + [fit.best_errors[ + fit.best_parameters()\ + .index(f'peak{i+1}_center')] + for i in range(num_peak)]) + unconstrained_fit_amplitudes = [ + fit.best_values[ + fit.best_parameters().index(f'peak{i+1}_amplitude')] + for i in range(num_peak)] + unconstrained_fit_amplitudes_errors = [ + fit.best_errors[ + fit.best_parameters().index(f'peak{i+1}_amplitude')] + for i in range(num_peak)] + unconstrained_fit_sigmas = [ + fit.best_values[ + fit.best_parameters().index(f'peak{i+1}_sigma')] + for i in range(num_peak)] + unconstrained_fit_sigmas_errors = [ + fit.best_errors[ + fit.best_parameters().index(f'peak{i+1}_sigma')] + for i in range(num_peak)] + unconstrained_best_fit = fit.best_fit + unconstrained_residuals = fit.residual + unconstrained_redchi = fit.redchi + unconstrained_success = fit.success + + return ( + uniform_fit_centers, uniform_fit_centers_errors, + uniform_fit_amplitudes, uniform_fit_amplitudes_errors, + uniform_fit_sigmas, uniform_fit_sigmas_errors, + uniform_best_fit, uniform_residuals, uniform_redchi, uniform_success, + unconstrained_fit_centers, unconstrained_fit_centers_errors, + unconstrained_fit_amplitudes, unconstrained_fit_amplitudes_errors, + unconstrained_fit_sigmas, unconstrained_fit_sigmas_errors, + unconstrained_best_fit, unconstrained_residuals, + unconstrained_redchi, unconstrained_success + ) diff --git a/CHAP/utils/scanparsers.py b/CHAP/utils/scanparsers.py index 307030b..bbd54c8 100755 --- a/CHAP/utils/scanparsers.py +++ b/CHAP/utils/scanparsers.py @@ -1114,8 +1114,6 @@ def __init__(self, spec_file_name, scan_number, detector_data_format=None): 'Unrecognized value for detector_data_format: ' + f'{detector_data_format}. Allowed values are: ' + ', '.join(self.detector_data_formats)) - print() - def init_detector_data_format(self): """Determine and set a value for the instance variable