Coverage for gpaw/utilities/folder.py: 93%

1import numpy as np 

2 

3from gpaw.gauss import Gauss 

4 

5 

6def x_y_xl(x, y, xl=None): 

7 X = np.array(x) 

8 assert len(X.shape) == 1 

9 Y = np.array(y) 

10 assert X.shape[0] == Y.shape[0] 

11 

12 if xl is None: 

13 Xl = np.unique(X) 

14 else: 

15 Xl = np.array(xl) 

16 assert len(Xl.shape) == 1 

17 return X, Y, Xl 

18 

19 

20class Lorentz: 

21 """Normalized Lorentzian distribution""" 

22 def __init__(self, width=0.08): 

23 self.dtype = float 

24 self.set_width(width) 

25 

26 def get(self, x, x0=0): 

27 return self.norm / ((x - x0)**2 + self.width2) 

28 

29 def set_width(self, width=0.08): 

30 self.norm = width / np.pi 

31 self.width2 = width**2 

32 self._fwhm = 2. * width 

33 

34 @property 

35 def fwhm(self): 

36 return self._fwhm 

37 

38 @fwhm.setter 

39 def fwhm(self, value): 

40 self.set_width(value / 2.) 

41 

42 

43class Voigt: 

44 """Voigt profile. 

45 

46 See https://reference.wolfram.com/language/ref/VoigtDistribution.html 

47 """ 

48 def __init__(self, width=0.08): 

49 self.dtype = float 

50 self.set_width(width) 

51 

52 def get(self, x, x0=0): 

53 from scipy.special import erfc 

54 argm = (-1j * (x - x0) + self.delta) * self.argpre 

55 argp = (1j * (x - x0) + self.delta) * self.argpre 

56 res = np.exp(argm**2) * erfc(argm) 

57 res += np.exp(argp**2) * erfc(argp) 

58 return res.real * self.prefactor 

59 

60 def set_width(self, width=0.08): 

61 """Width is interpreted as [delta, sigma]""" 

62 if not hasattr(width, '__iter__'): 

63 width = [width, width] 

64 self.delta, self.sigma = width 

65 self.width = np.linalg.norm(np.array(width)) 

66 self.prefactor = 1. / 2 / np.sqrt(2 * np.pi) / self.sigma 

67 self.argpre = 1. / np.sqrt(2) / self.sigma 

68 

69 """Full width at half maximum approximation after 

70 Olivero, J. J.; R. L. Longbothum 

71 "Empirical fits to the Voigt line width: A brief review" 

72 J. Quantitative Spectroscopy and Radiative Transfer. 17 (1977) 233-236 

73 """ 

74 fl = 2 * self.delta 

75 fg = np.sqrt(8 * np.log(2)) * self.sigma 

76 self._fwhm = 0.5346 * fl + np.sqrt(0.2166 * fl**2 + fg**2) 

77 

78 @property 

79 def fwhm(self): 

80 return self._fwhm 

81 

82 @fwhm.setter 

83 def fwhm(self, value): 

84 "setting fwhm is not uniquely defined" 

85 raise NotImplementedError 

86 

87 

88class ComplexLorentz: 

89 def __init__(self, width=0.08): 

90 self.dtype = complex 

91 self.set_width(width) 

92 

93 def get(self, x, x0): 

94 return 0.5 / x0 * (((x + x0) / ((x + x0)**2 + self.width2) - 

95 (x - x0) / ((x - x0)**2 + self.width2)) - 

96 1.0j * self.width * 

97 (1 / ((x + x0)**2 + self.width2) - 

98 1 / ((x - x0)**2 + self.width2)) 

99 ) 

100 

101 def set_width(self, width=0.08): 

102 self.width = width 

103 self.width2 = width**2 

104 

105 

106class ComplexGauss: 

107 def __init__(self, width=0.08): 

108 self.dtype = complex 

109 self.set_width(width) 

110 

111 def get(self, x, x0): 

112 from scipy.special import dawsn 

113 return 0.5 / x0 * (2 * self.wm1 * 

114 (dawsn((x + x0) * self.wm1) - 

115 dawsn((x - x0) * self.wm1)) - 

116 1.0j * self.norm * 

117 (np.exp(-((x + x0) * self.wm1)**2) - 

118 np.exp(-((x - x0) * self.wm1)**2)) 

119 ) 

120 

121 def set_width(self, width=0.08): 

122 self.norm = 1. / width * np.sqrt(np.pi / 2) 

123 self.wm1 = np.sqrt(.5) / width 

124 

125 

126class LorentzPole: 

127 """Pole of a damped harmonic oscillator. 

128 

129 See: e.g. Frederick Wooten, Optical properties of solids, 

130 Academic Press Inc. (1972)""" 

131 def __init__(self, width=0.08, imag=True): 

132 self.dtype = float 

133 if imag: 

134 self.get = self.get_imaginary 

135 else: 

136 self.get = self.get_real 

137 self.set_width(width) 

138 

139 def get_real(self, x, x0): 

140 return (x0**2 - x**2) / ( 

141 (x0**2 - x**2)**2 + self.width**2 * x**2) 

142 

143 def get_imaginary(self, x, x0): 

144 return self.width * x / ( 

145 (x0**2 - x**2)**2 + self.width**2 * x**2) 

146 

147 def set_width(self, width=0.08): 

148 self.width = width 

149 

150 

151class Folder: 

152 """Fold a function with normalised Gaussians or Lorentzians 

153 

154 Example: fold a function y(x) by Lorentzians of 0.2 width 

155 

156 >>> xlist = np.linspace(-1.0, 1.0, 101) # list of x-values 

157 >>> ylist = np.zeros_like(xlist) # corresponding list of y(x)-values 

158 >>> ylist[50] = 1.0 

159 >>> from gpaw.utilities.folder import Folder 

160 >>> fxlist, fylist = Folder(width=0.2, folding='Lorentz').fold( 

161 ... xlist, ylist, dx=0.1, xmin=-1, xmax=1) 

162 """ 

163 def __init__(self, width, 

164 folding='Gauss'): 

165 self.width = width 

166 if folding == 'Gauss': 

167 self.func = Gauss(width) 

168 elif folding == 'Lorentz': 

169 self.func = Lorentz(width) 

170 elif folding == 'ComplexGauss': 

171 self.func = ComplexGauss(width) 

172 elif folding == 'ComplexLorentz': 

173 self.func = ComplexLorentz(width) 

174 elif folding == 'RealLorentzPole': 

175 self.func = LorentzPole(width, imag=False) 

176 elif folding == 'ImaginaryLorentzPole': 

177 self.func = LorentzPole(width, imag=True) 

178 elif folding == 'Voigt': 

179 self.func = Voigt(width) 

180 else: 

181 raise RuntimeError('unknown folding "' + folding + '"') 

182 

183 def x_lim(self, x, dx=None, xmin=None, xmax=None): 

184 try: 

185 w = self.func.width 

186 except AttributeError: 

187 w = self.width 

188 

189 if xmin is None: 

190 xmin = np.min(x) - 4 * w 

191 if xmax is None: 

192 xmax = np.max(x) + 4 * w 

193 if dx is None: 

194 try: 

195 dx = self.func.width / 4. 

196 except AttributeError: 

197 dx = self.width / 4. 

198 

199 return dx, xmin, xmax 

200 

201 def fold(self, x, y, dx=None, xmin=None, xmax=None): 

202 X = np.array(x) 

203 assert len(X.shape) == 1 

204 Y = np.array(y) 

205 assert X.shape[0] == Y.shape[0] 

206 

207 dx, xmin, xmax = self.x_lim(X, dx, xmin, xmax) 

208 

209 xl = np.arange(xmin, xmax + 0.5 * dx, dx) 

210 

211 return self.fold_values(x, y, xl) 

212 

213 def fold_values(self, x, y, xl=None): 

214 X, Y, Xl = x_y_xl(x, y, xl) 

215 

216 # weight matrix 

217 weightm = np.empty((Xl.shape[0], X.shape[0]), 

218 dtype=self.func.dtype) 

219 for i, x in enumerate(X): 

220 weightm[:, i] = self.func.get(Xl, x) 

221 

222 yl = np.tensordot(weightm, Y, axes=(1, 0)) 

223 

224 return Xl, yl 

225 

226 def varing_fold(self, x, y, width2, x1, x2, 

227 dx=None, xmin=None, xmax=None): 

228 

229 X = np.array(x) 

230 assert len(X.shape) == 1 

231 Y = np.array(y) 

232 assert X.shape[0] == Y.shape[0] 

233 

234 dx, xmin, xmax = self.x_lim(X, dx, xmin, xmax) 

235 

236 xl = np.arange(xmin, xmax + 0.5 * dx, dx) 

237 return self.varing_fold_values(x, y, xl, width2, x1, x2) 

238 

239 def varing_fold_values(self, x, y, xl, width2, x1, x2,): 

240 

241 X, Y, Xl = x_y_xl(x, y, xl) 

242 

243 width1 = self.width 

244 

245 weightm = np.empty((Xl.shape[0], X.shape[0]), 

246 dtype=self.func.dtype) 

247 

248 for i, x in enumerate(X): 

249 if x < x1: 

250 width_eff = width1 

251 elif x < x2: 

252 width_eff = (width1 + (x - x1) * 

253 (width2 - width1) / (x2 - x1)) 

254 elif x >= x2: 

255 width_eff = width2 

256 self.func.set_width(width_eff) 

257 

258 weightm[:, i] = self.func.get(Xl, x) 

259 

260 yl = np.tensordot(weightm, Y, axes=(1, 0)) 

261 

262 return Xl, yl