Coverage for larch/xafs/mback.py: 8%
179 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1#!/usr/bin/env python
2"""
3 XAFS MBACK normalization algorithms.
4"""
5import numpy as np
6from scipy.special import erfc
8from xraydb import (xray_edge, xray_line, xray_lines,
9 f1_chantler, f2_chantler, guess_edge,
10 atomic_number, atomic_symbol)
11from lmfit import Parameter, Parameters, minimize
13from larch import Group, isgroup, parse_group_args, Make_CallArgs
15from larch.math import index_of, index_nearest, remove_dups, remove_nans2
17from .xafsutils import set_xafsGroup, TINY_ENERGY
18from .pre_edge import find_e0, preedge, pre_edge
20MAXORDER = 6
22def find_xray_line(z, edge):
23 """
24 Finds most intense X-ray emission line energy for a given element and edge.
25 """
26 intensity = 0
27 line = ''
28 for key, value in xray_lines(z).items() :
29 if value.initial_level == edge.upper():
30 if value.intensity > intensity:
31 intensity = value.intensity
32 line = key
33 return xray_line(z, line[:-1])
35def match_f2(p, en=0, mu=1, f2=1, e0=0, em=0, weight=1, theta=1, order=None,
36 leexiang=False):
37 """
38 Objective function for matching mu(E) data to tabulated f''(E) using the MBACK
39 algorithm and, optionally, the Lee & Xiang extension.
40 """
41 pars = p.valuesdict()
42 eoff = en - e0
44 norm = p['a']*erfc((en-em)/p['xi']) + p['c0'] # erfc function + constant term
45 for i in range(order): # successive orders of polynomial
46 attr = 'c%d' % (i + 1)
47 if attr in p:
48 norm += p[attr] * eoff**(i + 1)
49 func = (f2 + norm - p['s']*mu) * theta / weight
50 if leexiang:
51 func = func / p['s']*mu
52 return func
55def mback(energy, mu=None, group=None, z=None, edge='K', e0=None, pre1=None, pre2=-50,
56 norm1=100, norm2=None, order=3, leexiang=False, tables='chantler', fit_erfc=False,
57 return_f1=False, _larch=None):
58 """
59 Match mu(E) data for tabulated f''(E) using the MBACK algorithm and,
60 optionally, the Lee & Xiang extension
62 Arguments
63 ----------
64 energy: array of x-ray energies, in eV.
65 mu: array of mu(E).
66 group: output group.
67 z: atomic number of the absorber.
68 edge: x-ray absorption edge (default 'K')
69 e0: edge energy, in eV. If None, it will be determined here.
70 pre1: low E range (relative to e0) for pre-edge region.
71 pre2: high E range (relative to e0) for pre-edge region.
72 norm1: low E range (relative to e0) for post-edge region.
73 norm2: high E range (relative to e0) for post-edge region.
74 order: order of the legendre polynomial for normalization.
75 (default=3, min=0, max=5).
76 leexiang: boolean (default False) to use the Lee & Xiang extension.
77 tables: tabulated scattering factors: 'chantler' [deprecated]
78 fit_erfc: boolean (default False) to fit parameters of error function.
79 return_f1: boolean (default False) to include the f1 array in the group.
82 Returns
83 -------
84 None
86 The following attributes will be written to the output group:
87 group.f2: tabulated f2(E).
88 group.f1: tabulated f1(E) (if 'return_f1' is True).
89 group.fpp: mback atched spectrum.
90 group.edge_step: edge step of spectrum.
91 group.norm: normalized spectrum.
92 group.mback_params: group of parameters for the minimization.
94 Notes:
95 Chantler tables is now used, with Cromer-Liberman no longer supported.
96 References:
97 * MBACK (Weng, Waldo, Penner-Hahn): http://dx.doi.org/10.1086/303711
98 * Lee and Xiang: http://dx.doi.org/10.1088/0004-637X/702/2/970
99 * Cromer-Liberman: http://dx.doi.org/10.1063/1.1674266
100 * Chantler: http://dx.doi.org/10.1063/1.555974
101 """
102 order = max(min(order, MAXORDER), 0)
104 ### implement the First Argument Group convention
105 energy, mu, group = parse_group_args(energy, members=('energy', 'mu'),
106 defaults=(mu,), group=group,
107 fcn_name='mback')
108 if len(energy.shape) > 1:
109 energy = energy.squeeze()
110 if len(mu.shape) > 1:
111 mu = mu.squeeze()
113 if _larch is not None:
114 group = set_xafsGroup(group, _larch=_larch)
116 energy = remove_dups(energy, tiny=TINY_ENERGY)
117 if energy.size <= 1:
118 raise ValueError("energy array must have at least 2 points")
119 if e0 is None or e0 < energy[1] or e0 > energy[-2]:
120 e0 = find_e0(energy, mu, group=group)
122 ie0 = index_nearest(energy, e0)
123 e0 = energy[ie0]
125 pre1_input = pre1
126 norm2_input = norm2
128 if pre1 is None: pre1 = min(energy) - e0
129 if norm2 is None: norm2 = max(energy) - e0
130 if norm2 < 0: norm2 = max(energy) - e0 - norm2
131 pre1 = max(pre1, (min(energy) - e0))
132 norm2 = min(norm2, (max(energy) - e0))
134 if pre1 > pre2:
135 pre1, pre2 = pre2, pre1
136 if norm1 > norm2:
137 norm1, norm2 = norm2, norm1
139 p1 = index_of(energy, pre1+e0)
140 p2 = index_nearest(energy, pre2+e0)
141 n1 = index_nearest(energy, norm1+e0)
142 n2 = index_of(energy, norm2+e0)
143 if p2 - p1 < 2:
144 p2 = min(len(energy), p1 + 2)
145 if n2 - n1 < 2:
146 p2 = min(len(energy), p1 + 2)
148 ## theta is a boolean array indicating the
149 ## energy values considered for the fit.
150 ## theta=1 for included values, theta=0 for excluded values.
151 theta = np.zeros_like(energy, dtype='int')
152 theta[p1:(p2+1)] = 1
153 theta[n1:(n2+1)] = 1
155 ## weights for the pre- and post-edge regions, as defined in the MBACK paper (?)
156 weight = np.ones_like(energy, dtype=float)
157 weight[p1:(p2+1)] = np.sqrt(np.sum(weight[p1:(p2+1)]))
158 weight[n1:(n2+1)] = np.sqrt(np.sum(weight[n1:(n2+1)]))
160 ## get the f'' function from CL or Chantler
161 f1 = f1_chantler(z, energy)
162 f2 = f2_chantler(z, energy)
163 group.f2 = f2
164 if return_f1:
165 group.f1 = f1
167 em = find_xray_line(z, edge).energy # erfc centroid
169 params = Parameters()
170 params.add(name='s', value=1.0, vary=True) # scale of data
171 params.add(name='xi', value=50.0, vary=False, min=0) # width of erfc
172 params.add(name='a', value=0.0, vary=False) # amplitude of erfc
173 if fit_erfc:
174 params['a'].vary = True
175 params['a'].value = 0.5
176 params['xi'].vary = True
178 for i in range(order+1): # polynomial coefficients
179 params.add(name='c%d' % i, value=0, vary=True)
181 out = minimize(match_f2, params, method='leastsq',
182 gtol=1.e-5, ftol=1.e-5, xtol=1.e-5, epsfcn=1.e-5,
183 kws = dict(en=energy, mu=mu, f2=f2, e0=e0, em=em,
184 order=order, weight=weight, theta=theta, leexiang=leexiang))
186 opars = out.params.valuesdict()
187 eoff = energy - e0
189 norm_function = opars['a']*erfc((energy-em)/opars['xi']) + opars['c0']
190 for i in range(order):
191 attr = 'c%d' % (i + 1)
192 if attr in opars:
193 norm_function += opars[attr]* eoff**(i + 1)
195 group.e0 = e0
196 group.fpp = opars['s']*mu - norm_function
197 # calculate edge step and normalization from f2 + norm_function
198 pre_f2 = preedge(energy, group.f2+norm_function, e0=e0, pre1=pre1,
199 pre2=pre2, norm1=norm1, norm2=norm2, nnorm=2, nvict=0)
200 group.edge_step = pre_f2['edge_step'] / opars['s']
201 group.norm = (opars['s']*mu - pre_f2['pre_edge']) / pre_f2['edge_step']
202 group.mback_details = Group(params=opars, pre_f2=pre_f2,
203 f2_scaled=opars['s']*f2,
204 norm_function=norm_function)
208def f2norm(params, en=1, mu=1, f2=1, weights=1):
210 """
211 Objective function for matching mu(E) data to tabulated f''(E)
212 """
213 p = params.valuesdict()
214 model = (p['offset'] + p['slope']*en + f2) * p['scale']
215 return weights * (model - mu)
217@Make_CallArgs(["energy","mu"])
218def mback_norm(energy, mu=None, group=None, z=None, edge='K', e0=None,
219 pre1=None, pre2=None, norm1=None, norm2=None, nnorm=None, nvict=1,
220 _larch=None):
221 """
222 simplified version of MBACK to Match mu(E) data for tabulated f''(E)
223 for normalization
225 Arguments:
226 energy, mu: arrays of energy and mu(E)
227 group: output group (and input group for e0)
228 z: Z number of absorber
229 e0: edge energy
230 pre1: low E range (relative to E0) for pre-edge fit
231 pre2: high E range (relative to E0) for pre-edge fit
232 norm1: low E range (relative to E0) for post-edge fit
233 norm2: high E range (relative to E0) for post-edge fit
234 nnorm: degree of polynomial (ie, nnorm+1 coefficients will be
235 found) for post-edge normalization curve fit to the
236 scaled f2. Default=1 (linear)
238 Returns:
239 group.norm_poly: normalized mu(E) from pre_edge()
240 group.norm: normalized mu(E) from this method
241 group.mback_mu: tabulated f2 scaled and pre_edge added to match mu(E)
242 group.mback_params: Group of parameters for the minimization
244 References:
245 * MBACK (Weng, Waldo, Penner-Hahn): http://dx.doi.org/10.1086/303711
246 * Chantler: http://dx.doi.org/10.1063/1.555974
247 """
248 ### implement the First Argument Group convention
249 energy, mu, group = parse_group_args(energy, members=('energy', 'mu'),
250 defaults=(mu,), group=group,
251 fcn_name='mback')
252 if len(energy.shape) > 1:
253 energy = energy.squeeze()
254 if len(mu.shape) > 1:
255 mu = mu.squeeze()
257 if _larch is not None:
258 group = set_xafsGroup(group, _larch=_larch)
259 group.norm_poly = group.norm*1.0
261 if z is not None: # need to run find_e0:
262 e0_nominal = xray_edge(z, edge).energy
263 if e0 is None:
264 e0 = getattr(group, 'e0', None)
265 if e0 is None:
266 find_e0(energy, mu, group=group)
267 e0 = group.e0
269 atsym = None
270 if z is None or z < 2:
271 atsym, edge = guess_edge(group.e0)
272 z = atomic_number(atsym)
273 if atsym is None and z is not None:
274 atsym = atomic_symbol(z)
276 if getattr(group, 'pre_edge_details', None) is None: # pre_edge never run
277 pre_edge(group, pre1=pre1, pre2=pre2, nvict=nvict,
278 norm1=norm1, norm2=norm2, e0=e0, nnorm=nnorm)
279 if pre1 is None:
280 pre1 = group.pre_edge_details.pre1
281 if pre2 is None:
282 pre2 = group.pre_edge_details.pre2
283 if nvict is None:
284 nvict = group.pre_edge_details.nvict
285 if norm1 is None:
286 norm1 = group.pre_edge_details.norm1
287 if norm2 is None:
288 norm2 = group.pre_edge_details.norm2
289 if nnorm is None:
290 nnorm = group.pre_edge_details.nnorm
292 mu_pre = mu - group.pre_edge
293 f2 = f2_chantler(z, energy)
295 weights = np.ones(len(energy))*1.0
297 if norm2 is None:
298 norm2 = max(energy) - e0
299 if norm2 < 0:
300 norm2 = (max(energy) - e0) - norm2
302 # avoid l2 and higher edges
303 if edge.lower().startswith('l'):
304 if edge.lower() == 'l3':
305 e_l2 = xray_edge(z, 'L2').energy
306 norm2 = min(norm2, e_l2-e0)
307 elif edge.lower() == 'l2':
308 e_l1 = xray_edge(z, 'L1').energy
309 norm2 = min(norm2, e_l1-e0)
311 ipre2 = index_of(energy, e0+pre2)
312 inor1 = index_of(energy, e0+norm1)
313 inor2 = index_of(energy, e0+norm2) + 1
317 weights[ipre2:] = 0.0
318 weights[inor1:inor2] = np.linspace(0.1, 1.0, inor2-inor1)
320 params = Parameters()
321 params.add(name='slope', value=0.0, vary=True)
322 params.add(name='offset', value=-f2[0], vary=True)
323 params.add(name='scale', value=f2[-1], vary=True)
325 out = minimize(f2norm, params, method='leastsq',
326 gtol=1.e-5, ftol=1.e-5, xtol=1.e-5, epsfcn=1.e-5,
327 kws = dict(en=energy, mu=mu_pre, f2=f2, weights=weights))
329 p = out.params.valuesdict()
331 model = (p['offset'] + p['slope']*energy + f2) * p['scale']
333 group.mback_mu = model + group.pre_edge
335 pre_f2 = preedge(energy, model, nnorm=nnorm, nvict=nvict, e0=e0,
336 pre1=pre1, pre2=pre2, norm1=norm1, norm2=norm2)
338 step_new = pre_f2['edge_step']
340 group.edge_step_poly = group.edge_step
341 group.edge_step_mback = step_new
342 group.norm_mback = mu_pre / step_new
345 group.mback_params = Group(e0=e0, pre1=pre1, pre2=pre2, norm1=norm1,
346 norm2=norm2, nnorm=nnorm, fit_params=p,
347 fit_weights=weights, model=model, f2=f2,
348 pre_f2=pre_f2, atsym=atsym, edge=edge)
350 if (abs(step_new - group.edge_step)/(1.e-13+group.edge_step)) > 0.75:
351 print("Warning: mback edge step failed....")
352 else:
353 group.edge_step = step_new
354 group.norm = group.norm_mback