Coverage for larch/xafs/diffkk.py: 15%

136 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-10-16 21:04 +0000

1import time 

2import numpy as np 

3from scipy.special import erfc 

4 

5from larch import Group 

6from larch.math import interp 

7from larch.utils.show import show 

8 

9from .mback import mback 

10 

11pi = np.pi 

12 

13TINY = 1e-20 

14FOPI = 4/pi 

15MAXORDER = 6 

16 

17## 

18## to test the advantage of the vectorized MacLaurin series algorithm: 

19## 

20# import time 

21# a=read_ascii('/home/bruce/git/demeter/examples/cu/cu10k.dat') 

22# 

23# b=diffkk(a.energy, a.xmu, e0=8979, z=29, order=4, form='mback') 

24# start = time.clock() 

25# for i in range(10): 

26# b.kktrans() 

27# endfor 

28# finish = time.clock() 

29# print(finish - start) 

30# 

31# start = time.clock() 

32# for i in range(10): 

33# b.kktrans(how='sca') 

34# endfor 

35# finish = time.clock() 

36# print(finish - start) 

37## 

38## I got 7.6 seconds and 76.1 seconds for 10 iterations of each. Awesome! 

39## 

40 

41 

42### 

43### These are the scalar forms of the MacLaurin series algorithm as originally coded up by Matt 

44### see https://github.com/newville/ifeffit/blob/master/src/diffkk/kkmclr.f, 

45### https://github.com/newville/ifeffit/blob/master/src/diffkk/kkmclf.f, 

46### and https://gist.github.com/maurov/33997083e96ab4036fe7 

47### The last one is a direct translation of Fortan to python, written by Matt, and used in 

48### CARD (http://www.esrf.eu/computing/scientific/CARD/CARD.html) 

49### 

50### See below for vector forms (much faster!) 

51### 

52def kkmclf_sca(e, finp): 

53 """ 

54 forward (f'->f'') kk transform, using mclaurin series algorithm 

55 

56 arguments: 

57 npts size of arrays to consider 

58 e energy array *must be on an even grid with an even number of points* [npts] (in) 

59 finp f' array [npts] (in) 

60 fout f'' array [npts] (out) 

61 notes fopi = 4/pi 

62 """ 

63 npts = len(e) 

64 if npts != len(finp): 

65 Exception("Input arrays not of same length in kkmclf_sca") 

66 if npts < 2: 

67 Exception("Array too short in kkmclf_sca") 

68 if npts % 2: 

69 Exception("Array has an odd number of elements in kkmclf_sca") 

70 fout = [0.0]*npts 

71 if npts >= 2: 

72 factor = FOPI * (e[npts-1] - e[0]) / (npts - 1) 

73 nptsk = npts / 2 

74 for i in range(npts): 

75 fout[i] = 0.0 

76 ei2 = e[i]*e[i] 

77 ioff = i%2 - 1 

78 for k in range(nptsk): 

79 j = k + k + ioff 

80 de2 = e[j]*e[j] - ei2 

81 if abs(de2) <= TINY: 

82 de2 = TINY 

83 fout[i] = fout[i] + finp[j]/de2 

84 fout[i] *= factor*e[i] 

85 return fout 

86 

87def kkmclr_sca(e, finp): 

88 """ 

89 reverse (f''->f') kk transform, using maclaurin series algorithm 

90 

91 arguments: 

92 npts size of arrays to consider 

93 e energy array *must be on a even grid with an even number of points* [npts] (in) 

94 finp f'' array [npts] (in) 

95 fout f' array [npts] (out) 

96 m newville jan 1997 

97 """ 

98 npts = len(e) 

99 if npts != len(finp): 

100 Exception("Input arrays not of same length in kkmclr") 

101 if npts < 2: 

102 Exception("Array too short in kkmclr") 

103 if npts % 2: 

104 Exception("Array has an odd number of elements in kkmclr") 

105 fout = [0.0]*npts 

106 

107 factor = -FOPI * (e[npts-1] - e[0]) / (npts - 1) 

108 nptsk = npts / 2 

109 for i in range(npts): 

110 fout[i] = 0.0 

111 ei2 = e[i]*e[i] 

112 ioff = i%2 - 1 

113 for k in range(nptsk): 

114 j = k + k + ioff 

115 de2 = e[j]*e[j] - ei2 

116 if abs(de2) <= TINY: 

117 de2 = TINY 

118 fout[i] = fout[i] + e[j]*finp[j]/de2 

119 fout[i] *= factor 

120 return fout 

121 

122### 

123### These are vector forms of the MacLaurin series algorithm, adapted from Matt's code by Bruce 

124### They are about an order of magnitude faster. 

125### 

126 

127def kkmclf(e, finp): 

128 """ 

129 forward (f'->f'') kk transform, using maclaurin series algorithm 

130 

131 arguments: 

132 e energy array *must be on an even grid with an even number of points* [npts] (in) 

133 finp f'' array [npts] (in) 

134 fout f' array [npts] (out) 

135 """ 

136 npts = len(e) 

137 if npts != len(finp): 

138 Exception("Input arrays not of same length for diff KK transform in kkmclr") 

139 if npts < 2: 

140 Exception("Array too short for diff KK transform in kkmclr") 

141 if npts % 2: 

142 Exception("Array has an odd number of elements for diff KK transform in kkmclr") 

143 

144 fout = np.zeros(npts) 

145 factor = FOPI * (e[-1] - e[0]) / (npts-1) 

146 ei2 = e**2 

147 ioff = np.mod(np.arange(npts), 2) - 1 

148 

149 nptsk = int(npts/2) 

150 k = np.arange(nptsk) 

151 

152 for i in range(npts): 

153 j = 2*k + ioff[i] 

154 de2 = e[j]**2 - ei2[i] 

155 fout[i] = sum(finp[j]/de2) 

156 

157 fout = fout * factor 

158 return fout 

159 

160 

161def kkmclr(e, finp): 

162 """ 

163 reverse (f''->f') kk transform, using maclaurin series algorithm 

164 

165 arguments: 

166 e energy array *must be on an even grid with an even number of points* [npts] (in) 

167 finp f' array [npts] (in) 

168 fout f'' array [npts] (out) 

169 """ 

170 npts = len(e) 

171 if npts != len(finp): 

172 Exception("Input arrays not of same for length diff KK transform in kkmclr") 

173 if npts < 2: 

174 Exception("Array too short for diff KK transform in kkmclr") 

175 if npts % 2: 

176 Exception("Array has an odd number of elements for diff KK transform in kkmclr") 

177 

178 fout = np.zeros(npts) 

179 factor = -FOPI * (e[-1] - e[0]) / (npts-1) 

180 ei2 = e**2 

181 ioff = np.mod(np.arange(npts), 2) - 1 

182 

183 nptsk = int(npts/2) 

184 k = np.arange(nptsk) 

185 

186 for i in range(npts): 

187 j = 2*k + ioff[i] 

188 de2 = e[j]**2 - ei2[i] 

189 fout[i] = sum(e[j]*finp[j]/de2) 

190 

191 fout = fout * factor 

192 return fout 

193 

194 

195class diffKKGroup(Group): 

196 """ 

197 A Larch Group for generating f'(E) and f"(E) from a XAS measurement of mu(E). 

198 """ 

199 

200 def __init__(self, energy=None, mu=None, z=None, edge='K', mback_kws=None, **kws): 

201 kwargs = dict(name='diffKK') 

202 kwargs.update(kws) 

203 Group.__init__(self, **kwargs) 

204 self.energy = energy 

205 self.mu = mu 

206 self.z = z 

207 self.edge = edge 

208 self.mback_kws = mback_kws 

209 def __repr__(self): 

210 return '<diffKK Group>' 

211 

212 

213# e0=None, z=None, edge=None, order=3, form='mback', whiteline=False, how=None 

214 def kk(self, energy=None, mu=None, z=None, edge='K', how='scalar', mback_kws=None): 

215 """ 

216 Convert mu(E) data into f'(E) and f"(E). f"(E) is made by 

217 matching mu(E) to the tabulated values of the imaginary part 

218 of the scattering factor (Cromer-Liberman), f'(E) is then 

219 obtained by performing a differential Kramers-Kronig transform 

220 on the matched f"(E). 

221 

222 Attributes 

223 energy: energy array 

224 mu: array with mu(E) data 

225 z: Z number of absorber 

226 edge: absorption edge, usually 'K' or 'L3' 

227 mback_kws: arguments for the mback algorithm 

228 

229 Returns 

230 self.f1, self.f2: CL values over on the input energy grid 

231 self.fp, self.fpp: matched and KK transformed data on the input energy grid 

232 

233 References: 

234 * Cromer-Liberman: http://dx.doi.org/10.1063/1.1674266 

235 * KK computation: Ohta and Ishida, Applied Spectroscopy 42:6 (1988) 952-957 

236 * diffKK implementation: http://dx.doi.org/10.1103/PhysRevB.58.11215 

237 * MBACK (Weng, Waldo, Penner-Hahn): http://dx.doi.org/10.1086/303711 

238 * Lee and Xiang: http://dx.doi.org/10.1088/0004-637X/702/2/970 

239 

240 """ 

241 

242 if type(energy).__name__ == 'ndarray': self.energy = energy 

243 if type(mu).__name__ == 'ndarray': self.mu = mu 

244 if z != None: self.z = z 

245 if edge != None: self.edge = edge 

246 if mback_kws != None: self.mback_kws = mback_kws 

247 

248 if self.z == None: 

249 Exception("Z for absorber not provided for diffKK") 

250 if self.edge == None: 

251 Exception("absorption edge not provided for diffKK") 

252 

253 mb_kws = dict(order=3, z=self.z, edge=self.edge, e0=None, 

254 leexiang=False, tables='chantler', 

255 fit_erfc=False, return_f1=True) 

256 if self.mback_kws is not None: 

257 mb_kws.update(self.mback_kws) 

258 

259 start = time.monotonic() 

260 

261 mback(self.energy, self.mu, group=self, **mb_kws) 

262 

263 ## interpolate matched data onto an even grid with an even number of elements (about 1 eV) 

264 npts = int(self.energy[-1] - self.energy[0]) + (int(self.energy[-1] - self.energy[0])%2) 

265 self.grid = np.linspace(self.energy[0], self.energy[-1], npts) 

266 fpp = interp(self.energy, self.f2-self.fpp, self.grid, fill_value=0.0) 

267 

268 ## do difference KK 

269 if repr(how).startswith('sca'): 

270 fp = kkmclr_sca(self.grid, fpp) 

271 else: 

272 fp = kkmclr(self.grid, fpp) 

273 

274 ## interpolate back to original grid and add diffKK result to f1 to make fp array 

275 self.fp = self.f1 + interp(self.grid, fp, self.energy, fill_value=0.0) 

276 

277 ## clean up group 

278 #for att in ('normalization_function', 'weight', 'grid'): 

279 # if hasattr(self, att): delattr(self, att) 

280 finish = time.monotonic() 

281 self.time_elapsed = float(finish-start) 

282 

283def diffkk(energy=None, mu=None, z=None, edge='K', mback_kws=None, **kws): 

284 """ 

285 Make a diffKK group given mu(E) data 

286 

287 Attributes 

288 energy: energy array 

289 mu: array with mu(E) data 

290 z: Z number of absorber 

291 edge: absorption edge, usually 'K' or 'L3' 

292 mback_kws: arguments for the mback algorithm 

293 """ 

294 return diffKKGroup(energy=energy, mu=mu, z=z, mback_kws=mback_kws)