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

1#!/usr/bin/env python 

2""" 

3 XAFS MBACK normalization algorithms. 

4""" 

5import numpy as np 

6from scipy.special import erfc 

7 

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 

12 

13from larch import Group, isgroup, parse_group_args, Make_CallArgs 

14 

15from larch.math import index_of, index_nearest, remove_dups, remove_nans2 

16 

17from .xafsutils import set_xafsGroup, TINY_ENERGY 

18from .pre_edge import find_e0, preedge, pre_edge 

19 

20MAXORDER = 6 

21 

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

34 

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 

43 

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 

53 

54 

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 

61 

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. 

80 

81 

82 Returns 

83 ------- 

84 None 

85 

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. 

93 

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) 

103 

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

112 

113 if _larch is not None: 

114 group = set_xafsGroup(group, _larch=_larch) 

115 

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) 

121 

122 ie0 = index_nearest(energy, e0) 

123 e0 = energy[ie0] 

124 

125 pre1_input = pre1 

126 norm2_input = norm2 

127 

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

133 

134 if pre1 > pre2: 

135 pre1, pre2 = pre2, pre1 

136 if norm1 > norm2: 

137 norm1, norm2 = norm2, norm1 

138 

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) 

147 

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 

154 

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

159 

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 

166 

167 em = find_xray_line(z, edge).energy # erfc centroid 

168 

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 

177 

178 for i in range(order+1): # polynomial coefficients 

179 params.add(name='c%d' % i, value=0, vary=True) 

180 

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

185 

186 opars = out.params.valuesdict() 

187 eoff = energy - e0 

188 

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) 

194 

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) 

205 

206 

207 

208def f2norm(params, en=1, mu=1, f2=1, weights=1): 

209 

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) 

216 

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 

224 

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) 

237 

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 

243 

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

256 

257 if _larch is not None: 

258 group = set_xafsGroup(group, _larch=_larch) 

259 group.norm_poly = group.norm*1.0 

260 

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 

268 

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) 

275 

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 

291 

292 mu_pre = mu - group.pre_edge 

293 f2 = f2_chantler(z, energy) 

294 

295 weights = np.ones(len(energy))*1.0 

296 

297 if norm2 is None: 

298 norm2 = max(energy) - e0 

299 if norm2 < 0: 

300 norm2 = (max(energy) - e0) - norm2 

301 

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) 

310 

311 ipre2 = index_of(energy, e0+pre2) 

312 inor1 = index_of(energy, e0+norm1) 

313 inor2 = index_of(energy, e0+norm2) + 1 

314 

315 

316 

317 weights[ipre2:] = 0.0 

318 weights[inor1:inor2] = np.linspace(0.1, 1.0, inor2-inor1) 

319 

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) 

324 

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

328 

329 p = out.params.valuesdict() 

330 

331 model = (p['offset'] + p['slope']*energy + f2) * p['scale'] 

332 

333 group.mback_mu = model + group.pre_edge 

334 

335 pre_f2 = preedge(energy, model, nnorm=nnorm, nvict=nvict, e0=e0, 

336 pre1=pre1, pre2=pre2, norm1=norm1, norm2=norm2) 

337 

338 step_new = pre_f2['edge_step'] 

339 

340 group.edge_step_poly = group.edge_step 

341 group.edge_step_mback = step_new 

342 group.norm_mback = mu_pre / step_new 

343 

344 

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) 

349 

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