From 6143f2d4ed77d8c94ff49e2fc51bbd560e2ce794 Mon Sep 17 00:00:00 2001 From: Rolf Verberg Date: Fri, 10 Nov 2023 12:42:13 -0500 Subject: [PATCH] fix: resolved EDD issues while analysing merrill-3229-f CHAP/edd/models.py: Fixed type error in background parameters CHAP/edd/processor.py: - Output DVL center in absolute coordinates - Display DVL relative to the correct center - Constrain peak centers to the actual energy range with some margin to allow for fitting without hitting boundary issues. - Continue unconstrained fit from uniform fit in Ceria calibration instead of starting with only the uniform peak centers CHAP/utils/fit.py: Add option to constrain peak centers to a given x-axis range in Fit:create_multipeak_model() --- CHAP/edd/models.py | 4 +- CHAP/edd/processor.py | 73 ++++++++++++++++----------- CHAP/utils/fit.py | 115 ++++++++++++++++++++++++++++-------------- 3 files changed, 123 insertions(+), 69 deletions(-) diff --git a/CHAP/edd/models.py b/CHAP/edd/models.py index a92e37a..377266f 100755 --- a/CHAP/edd/models.py +++ b/CHAP/edd/models.py @@ -405,7 +405,7 @@ class MCAElementCalibrationConfig(MCAElementConfig): tth_max: confloat(gt=0, allow_inf_nan=False) = 90.0 hkl_tth_tol: confloat(gt=0, allow_inf_nan=False) = 0.15 hkl_indices: Optional[conlist(item_type=conint(ge=0), min_items=1)] = [] - background: Optional[str] + background: Optional[Union[str, list]] tth_initial_guess: confloat(gt=0, le=tth_max, allow_inf_nan=False) = 5.0 slope_initial_guess: float = 1.0 intercept_initial_guess: float = 0.0 @@ -673,7 +673,7 @@ class MCAElementStrainAnalysisConfig(MCAElementConfig): tth_max: confloat(gt=0, allow_inf_nan=False) = 90.0 hkl_tth_tol: confloat(gt=0, allow_inf_nan=False) = 0.15 hkl_indices: Optional[conlist(item_type=conint(ge=0), min_items=1)] = [] - background: Optional[str] + background: Optional[Union[str, list]] peak_models: Union[ conlist(item_type=Literal['gaussian', 'lorentzian'], min_items=1), Literal['gaussian', 'lorentzian']] = 'gaussian' diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py index 74b0ca2..ad7473c 100755 --- a/CHAP/edd/processor.py +++ b/CHAP/edd/processor.py @@ -188,7 +188,7 @@ def measure_dvl(self, dvl = fit.best_values['sigma'] * detector.sigma_to_dvl_factor \ - dvl_config.sample_thickness detector.fit_amplitude = fit.best_values['amplitude'] - detector.fit_center = fit.best_values['center'] + detector.fit_center = scan_center + fit.best_values['center'] detector.fit_sigma = fit.best_values['sigma'] if detector.measurement_mode == 'manual': if interactive: @@ -226,7 +226,9 @@ def measure_dvl(self, ax.plot(x, masked_max, label='maximum (masked)') ax.plot(x, unmasked_sum, label='total (unmasked)') ax.axvspan( - -dvl / 2., dvl / 2., color='gray', alpha=0.5, + fit.best_values['center']- dvl/2., + fit.best_values['center'] + dvl/2., + color='gray', alpha=0.5, label=f'diffraction volume ({detector.measurement_mode})') ax.legend() plt.figtext( @@ -391,7 +393,7 @@ def calibrate(self, detector.hkl_indices = hkl_indices if save_figures: fig.savefig(os.path.join( - outputdir, + outputdir, f'{detector.detector_name}_calibration_fit_mask_hkls.png')) plt.close() self.logger.debug(f'tth_initial_guess = {detector.tth_initial_guess}') @@ -400,13 +402,13 @@ def calibrate(self, if detector.include_bin_ranges is None: raise ValueError( 'No value provided for include_bin_ranges. ' - 'Provide them in the MCA Ceria Calibration Configuration, ' + 'Provide them in the MCA Ceria Calibration Configuration ' 'or re-run the pipeline with the --interactive flag.') if detector.hkl_indices is None: raise ValueError( 'No value provided for hkl_indices. Provide them in ' - 'the detector\'s MCA Ceria Calibration Configuration, or' - ' re-run the pipeline with the --interactive flag.') + 'the detector\'s MCA Ceria Calibration Configuration or ' + 're-run the pipeline with the --interactive flag.') mca_mask = detector.mca_mask() fit_mca_energies = mca_bin_energies[mca_mask] fit_mca_intensities = mca_data[mca_mask] @@ -419,13 +421,18 @@ def calibrate(self, # Get the HKLs and lattice spacings that will be used for # fitting + # Restrict the range for the centers with some margin to + # prevent having centers near the edge of the fitting range + delta = 0.1 * (fit_mca_energies[-1]-fit_mca_energies[0]) + centers_range = ( + max(0.0, fit_mca_energies[0]-delta), fit_mca_energies[-1]+delta) fit_hkls = np.asarray([hkls[i] for i in detector.hkl_indices]) fit_ds = np.asarray([ds[i] for i in detector.hkl_indices]) c_1 = fit_hkls[:,0]**2 + fit_hkls[:,1]**2 + fit_hkls[:,2]**2 tth = detector.tth_initial_guess for iter_i in range(calibration_config.max_iter): - self.logger.debug(f'Tuning tth: iteration no. {iter_i}, ' - + f'starting tth value = {tth} ') + self.logger.debug( + f'Tuning tth: iteration no. {iter_i}, starting value = {tth} ') # Perform the uniform fit first @@ -434,36 +441,38 @@ def calibrate(self, fit_E0 = get_peak_locations(fit_ds, tth) # Run the uniform fit - uniform_fit = Fit(fit_mca_intensities, x=fit_mca_energies) - uniform_fit.create_multipeak_model( - fit_E0, fit_type='uniform', background=detector.background) - uniform_fit.fit() + fit = Fit(fit_mca_intensities, x=fit_mca_energies) + fit.create_multipeak_model( + fit_E0, fit_type='uniform', background=detector.background, + centers_range=centers_range) + fit.fit() # Extract values of interest from the best values for the # uniform fit parameters + uniform_best_fit = fit.best_fit + uniform_residual = fit.residual uniform_fit_centers = [ - uniform_fit.best_values[f'peak{i+1}_center'] + fit.best_values[f'peak{i+1}_center'] for i in range(len(fit_hkls))] - uniform_a = uniform_fit.best_values['scale_factor'] + uniform_a = fit.best_values['scale_factor'] uniform_strain = np.log( (uniform_a / calibration_config.material.lattice_parameters)) # CeO2 is cubic, so this is fine here. # Next, perform the unconstrained fit - # Use the peak locations found in the uniform fit as the + # Use the peak parameters from the uniform fit as the # initial guesses for peak locations in the unconstrained # fit - unconstrained_fit = Fit(fit_mca_intensities, x=fit_mca_energies) - unconstrained_fit.create_multipeak_model( - uniform_fit_centers, fit_type='unconstrained', - background=detector.background) - unconstrained_fit.fit() + fit.create_multipeak_model(fit_type='unconstrained') + fit.fit() # Extract values of interest from the best values for the # unconstrained fit parameters + unconstrained_best_fit = fit.best_fit + unconstrained_residual = fit.residual unconstrained_fit_centers = np.array( - [unconstrained_fit.best_values[f'peak{i+1}_center'] + [fit.best_values[f'peak{i+1}_center'] for i in range(len(fit_hkls))]) unconstrained_a = np.sqrt(c_1)*abs(get_peak_locations( unconstrained_fit_centers, tth)) @@ -505,13 +514,13 @@ def calibrate(self, axs[0,0].text(hkl_E, 1, str(fit_hkls[i])[1:-1], ha='right', va='top', rotation=90, transform=axs[0,0].get_xaxis_transform()) - axs[0,0].plot(fit_mca_energies, uniform_fit.best_fit, - label='Single Strain') - axs[0,0].plot(fit_mca_energies, unconstrained_fit.best_fit, - label='Unconstrained') + axs[0,0].plot(fit_mca_energies, uniform_best_fit, + label='Single Strain') + axs[0,0].plot(fit_mca_energies, unconstrained_best_fit, + label='Unconstrained') #axs[0,0].plot(fit_mca_energies, MISSING?, label='least squares') axs[0,0].plot(fit_mca_energies, fit_mca_intensities, - label='Flux-Corrected & Masked MCA Data') + label='Flux-Corrected & Masked MCA Data') axs[0,0].legend() # Lower left axes: fit residuals @@ -519,10 +528,10 @@ def calibrate(self, axs[1,0].set_xlabel('Energy (keV)') axs[1,0].set_ylabel('Residual (a.u)') axs[1,0].plot(fit_mca_energies, - uniform_fit.residual, + uniform_residual, label='Single Strain') axs[1,0].plot(fit_mca_energies, - unconstrained_fit.residual, + unconstrained_residual, label='Unconstrained') axs[1,0].legend() @@ -968,13 +977,17 @@ def linkdims(nxgroup, field_dims=[]): # 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) + fwhm_max=detector.fwhm_max, + centers_range=centers_range) fit.fit() uniform_fit_centers = [ fit.best_values[ @@ -1180,7 +1193,7 @@ def animate(i): path = os.path.join( outputdir, f'{detector.detector_name}_strainanalysis_' - 'unconstrained_fits.mp4') + 'unconstrained_fits.gif') ani.save(path) plt.close() diff --git a/CHAP/utils/fit.py b/CHAP/utils/fit.py index be84f55..8060279 100755 --- a/CHAP/utils/fit.py +++ b/CHAP/utils/fit.py @@ -862,9 +862,13 @@ def create_multipeak_model( # Third party modules from asteval import Interpreter + # Local modules + from CHAP.utils.general import is_num_pair + if centers_range is None: centers_range = (self._x[0], self._x[-1]) - elif not is_index_range(centers_range, ge=self._x[0], le=self._x[-1]): + elif (not is_num_pair(centers_range) or len(centers_range) != 2 + or centers_range[0] >= centers_range[1]): raise ValueError( f'Invalid parameter centers_range ({centers_range})') if self._model is not None: @@ -1021,9 +1025,8 @@ def create_multipeak_model( {'name': par_name, 'min': min_value}) elif len(index) == 1: parameter = parameters[index[0]] - _min = parameter.get('min', None) - if _min is None or _min < min_value: - parameter['min'] = min_value + _min = parameter.get('min', min_value) + parameter['min'] = max(_min, min_value) else: raise ValueError( 'Invalid parameters value in ' @@ -1045,14 +1048,47 @@ def create_multipeak_model( {'name': par_name, 'min': min_value}) elif len(index) == 1: parameter = parameters[index[0]] - _min = parameter.get('min', None) - if _min is None or _min < min_value: - parameter['min'] = min_value + _min = parameter.get('min', min_value) + parameter['min'] = max(_min, min_value) else: raise ValueError( 'Invalid parameters value in ' f'background model {name} ' f'({parameters})') + if name == 'gaussian': + if parameters is None: + parameters = { + 'name': 'center', + 'value': 0.5 * ( + centers_range[0] + centers_range[1]), + 'min': centers_range[0], + 'min': centers_range[1], + } + else: + index = [i for i, par in enumerate(parameters) + if par['name'] == 'center'] + if not len(index): + parameters.append({ + 'name': 'center', + 'value': 0.5 * ( + centers_range[0] + centers_range[1]), + 'min': centers_range[0], + 'max': centers_range[1], + }) + elif len(index) == 1: + parameter = parameters[index[0]] + if 'value' not in parameter: + parameter['value'] = 0.5 * ( + centers_range[0]+centers_range[1]) + _min = parameter.get('min', centers_range[0]) + parameter['min'] = max(_min, centers_range[0]) + _max = parameter.get('max', centers_range[1]) + parameter['max'] = min(_max, centers_range[1]) + else: + raise ValueError( + 'Invalid parameters value in ' + f'background model {name} ' + f'({parameters})') self.add_model( name, prefix=prefix, parameters=parameters, **model) @@ -1866,6 +1902,7 @@ def _renormalize(self): self._result.residual *= self._norm[1] def _reset_par_at_boundary(self): + fraction = 0.02 for name, par in self._parameters.items(): if par.vary: value = par.value @@ -1874,26 +1911,26 @@ def _reset_par_at_boundary(self): if np.isinf(_min): if not np.isinf(_max): if self._parameter_norms.get(name, False): - upp = _max-0.1*self._y_range + upp = _max - fraction*self._y_range elif _max == 0.0: - upp = _max-0.1 + upp = _max - fraction else: - upp = _max-0.1*abs(_max) + upp = _max - fraction*abs(_max) if value >= upp: par.set(value=upp) else: if np.isinf(_max): if self._parameter_norms.get(name, False): - low = _min + 0.1*self._y_range + low = _min + fraction*self._y_range elif _min == 0.0: - low = _min+0.1 + low = _min + fraction else: - low = _min + 0.1*abs(_min) + low = _min + fraction*abs(_min) if value <= low: par.set(value=low) else: - low = 0.9*_min + 0.1*_max - upp = 0.1*_min + 0.9*_max + low = (1.0-fraction)*_min + fraction*_max + upp = fraction*_min + (1.0-fraction)*_max if value <= low: par.set(value=low) if value >= upp: @@ -1932,6 +1969,7 @@ def __init__( # map dimensions if isinstance(ymap, (tuple, list, np.ndarray)): self._x = np.asarray(x) + ymap = np.asarray(ymap) elif HAVE_XARRAY and isinstance(ymap, xr.DataArray): if x is not None: logger.warning('Ignoring superfluous input x ({x})') @@ -2598,29 +2636,32 @@ def fit(self, **kwargs): except AttributeError: pass - if num_proc == 1: - # Perform the remaining fits serially - for n in range(1, self._map_dim): - self._fit(n, current_best_values, **kwargs) - else: - # Perform the remaining fits in parallel - num_fit = self._map_dim-1 - if num_proc > num_fit: - logger.warning( - f'The requested number of processors ({num_proc}) exceeds ' - f'the number of fits, num_proc reduced to {num_fit}') - num_proc = num_fit - num_fit_per_proc = 1 + if self._map_dim > 1: + if num_proc == 1: + # Perform the remaining fits serially + for n in range(1, self._map_dim): + self._fit(n, current_best_values, **kwargs) else: - num_fit_per_proc = round((num_fit)/num_proc) - if num_proc*num_fit_per_proc < num_fit: - num_fit_per_proc += 1 - num_fit_batch = min(num_fit_per_proc, 40) - with Parallel(n_jobs=num_proc) as parallel: - parallel( - delayed(self._fit_parallel) - (current_best_values, num_fit_batch, n_start, **kwargs) - for n_start in range(1, self._map_dim, num_fit_batch)) + # Perform the remaining fits in parallel + num_fit = self._map_dim-1 + if num_proc > num_fit: + logger.warning( + f'The requested number of processors ({num_proc}) ' + 'exceeds the number of fits, num_proc reduced to ' + f'{num_fit}') + num_proc = num_fit + num_fit_per_proc = 1 + else: + num_fit_per_proc = round((num_fit)/num_proc) + if num_proc*num_fit_per_proc < num_fit: + num_fit_per_proc += 1 + num_fit_batch = min(num_fit_per_proc, 40) + with Parallel(n_jobs=num_proc) as parallel: + parallel( + delayed(self._fit_parallel) + (current_best_values, num_fit_batch, n_start, + **kwargs) + for n_start in range(1, self._map_dim, num_fit_batch)) # Renormalize the initial parameters for external use if self._norm is not None and self._normalized: