Coverage for gpaw/response/hilbert.py: 100%

1import numpy as np 

2from gpaw.utilities.blas import mmm 

3 

4 

5class HilbertTransform: 

6 def __init__(self, omega_w, eta, timeordered=False, gw=False, 

7 blocksize=500): 

8 """Analytic Hilbert transformation using linear interpolation. 

9 

10 Hilbert transform:: 

11 

12 oo 

13 / 1 1 

14 |dw' (-------------- - --------------) S(w'). 

15 / w - w' + i eta w + w' + i eta 

16 0 

17 

18 With timeordered=True, you get:: 

19 

20 oo 

21 / 1 1 

22 |dw' (-------------- - --------------) S(w'). 

23 / w - w' - i eta w + w' + i eta 

24 0 

25 

26 With gw=True, you get:: 

27 

28 oo 

29 / 1 1 

30 |dw' (-------------- + --------------) S(w'). 

31 / w - w' + i eta w + w' + i eta 

32 0 

33 

34 """ 

35 

36 self.blocksize = blocksize 

37 

38 if timeordered: 

39 self.H_ww = self.H(omega_w, -eta) + self.H(omega_w, -eta, -1) 

40 elif gw: 

41 self.H_ww = self.H(omega_w, eta) - self.H(omega_w, -eta, -1) 

42 else: 

43 self.H_ww = self.H(omega_w, eta) + self.H(omega_w, -eta, -1) 

44 

45 def H(self, o_w, eta, sign=1): 

46 """Calculate transformation matrix. 

47 

48 With s=sign (+1 or -1):: 

49 

50 oo 

51 / dw' 

52 X (w, eta) = | ---------------- S(w'). 

53 s / s w - w' + i eta 

54 0 

55 

56 Returns H_ij so that X_i = np.dot(H_ij, S_j), where:: 

57 

58 X_i = X (omega_w[i]) and S_j = S(omega_w[j]) 

59 s 

60 """ 

61 

62 nw = len(o_w) 

63 H_ij = np.zeros((nw, nw), complex) 

64 do_j = o_w[1:] - o_w[:-1] 

65 for i, o in enumerate(o_w): 

66 d_j = o_w - o * sign 

67 y_j = 1j * np.arctan(d_j / eta) + 0.5 * np.log(d_j**2 + eta**2) 

68 y_j = (y_j[1:] - y_j[:-1]) / do_j 

69 H_ij[i, :-1] = 1 - (d_j[1:] - 1j * eta) * y_j 

70 H_ij[i, 1:] -= 1 - (d_j[:-1] - 1j * eta) * y_j 

71 return H_ij 

72 

73 def __call__(self, S_wx): 

74 """Inplace transform""" 

75 B_wx = S_wx.reshape((len(S_wx), -1)) 

76 nw, nx = B_wx.shape 

77 tmp_wx = np.zeros((nw, min(nx, self.blocksize)), complex) 

78 for x in range(0, nx, self.blocksize): 

79 b_wx = B_wx[:, x:x + self.blocksize] 

80 c_wx = tmp_wx[:, :b_wx.shape[1]] 

81 mmm(1.0, self.H_ww, 'N', b_wx, 'N', 0.0, c_wx) 

82 b_wx[:] = c_wx 

83 

84 

85class GWHilbertTransforms: 

86 """Helper class which wraps two transforms using contiguous array. 

87 

88 (This slightly speeds up things.)""" 

89 def __init__(self, omega_w, eta): 

90 self.htp = htp = HilbertTransform(omega_w, eta, gw=True) 

91 self.htm = htm = HilbertTransform(omega_w, -eta, gw=True) 

92 self._stacked_H_nww = np.array([htp.H_ww, htm.H_ww]) 

93 

94 def __call__(self, A_wGG): 

95 # (Note: This effectively duplicates the Hilbert call) 

96 nw = len(A_wGG) 

97 H_xw = self._stacked_H_nww.reshape(-1, nw) 

98 A_wy = A_wGG.reshape(nw, -1) 

99 tmp_xy = np.zeros((H_xw.shape[0], A_wy.shape[1]), complex) 

100 # gemm(1.0, A_wy, H_xw, 0.0, tmp_xy) 

101 mmm(1.0, H_xw, 'N', A_wy, 'N', 0.0, tmp_xy) 

102 return tmp_xy.reshape((2, *A_wGG.shape))