Skip to content

Commit

Permalink
fix: resolved EDD issues while analysing merrill-3229-f
Browse files Browse the repository at this point in the history 
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()
  • Loading branch information
@rolfverberg
rolfverberg committed Nov 10, 2023
1 parent e111758 commit 6143f2d
Show file tree
Hide file tree
Showing 3 changed files with 123 additions and 69 deletions.
4 changes: 2 additions & 2 deletions CHAP/edd/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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'
Expand Down
73 changes: 43 additions & 30 deletions CHAP/edd/processor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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}')
Expand All @@ -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]
Expand All @@ -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

Expand All @@ -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))
Expand Down Expand Up @@ -505,24 +514,24 @@ 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
axs[1,0].set_title('Fit Residuals')
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()

Expand Down Expand Up @@ -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[
Expand Down Expand Up @@ -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()

Expand Down
115 changes: 78 additions & 37 deletions CHAP/utils/fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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 '
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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})')
Expand Down Expand Up @@ -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:
Expand Down

0 comments on commit 6143f2d

Please sign in to comment.