Add CHAP.edd.utils.get_spectra_fits

Preparation for edd.LatticeParameterRefinementProcessor -- it will
share code with edd.StrainAnalysisProcessor.

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()

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
+    )

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