Coverage for pySDC/core/Lagrange.py: 100%

1import numpy as np 

2from scipy.special import roots_legendre 

3 

4 

5def computeFejerRule(n): 

6 """ 

7 Compute a Fejer rule of the first kind, using DFT (Waldvogel 2006) 

8 Inspired from quadpy (https://github.com/nschloe/quadpy @Nico_Schlömer) 

9 

10 Parameters 

11 ---------- 

12 n : int 

13 Number of points for the quadrature rule. 

14 

15 Returns 

16 ------- 

17 nodes : np.1darray(n) 

18 The nodes of the quadrature rule 

19 weights : np.1darray(n) 

20 The weights of the quadrature rule. 

21 """ 

22 # Initialize output variables 

23 n = int(n) 

24 nodes = np.empty(n, dtype=float) 

25 weights = np.empty(n, dtype=float) 

26 

27 # Compute nodes 

28 theta = np.arange(1, n + 1, dtype=float)[-1::-1] 

29 theta *= 2 

30 theta -= 1 

31 theta *= np.pi / (2 * n) 

32 np.cos(theta, out=nodes) 

33 

34 # Compute weights 

35 # -- Initial variables 

36 N = np.arange(1, n, 2) 

37 lN = len(N) 

38 m = n - lN 

39 K = np.arange(m) 

40 # -- Build v0 

41 v0 = np.concatenate([2 * np.exp(1j * np.pi * K / n) / (1 - 4 * K**2), np.zeros(lN + 1)]) 

42 # -- Build v1 from v0 

43 v1 = np.empty(len(v0) - 1, dtype=complex) 

44 np.conjugate(v0[:0:-1], out=v1) 

45 v1 += v0[:-1] 

46 # -- Compute inverse Fourier transform 

47 w = np.fft.ifft(v1) 

48 if max(w.imag) > 1.0e-15: 

49 raise ValueError(f'Max imaginary value to important for ifft: {max(w.imag)}') 

50 # -- Store weights 

51 weights[:] = w.real 

52 

53 return nodes, weights 

54 

55 

56class LagrangeApproximation(object): 

57 r""" 

58 Class approximating any function on a given set of points using barycentric 

59 Lagrange interpolation. 

60 

61 Let note :math:`(t_j)_{0\leq j<n}` the set of points, then any scalar 

62 function :math:`f` can be approximated by the barycentric formula : 

63 

64 .. math:: 

65 p(x) = 

66 \frac{\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}f_j} 

67 {\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}}, 

68 

69 where :math:`f_j=f(t_j)` and 

70 

71 .. math:: 

72 w_j = \frac{1}{\prod_{k\neq j}(x_j-x_k)} 

73 

74 are the barycentric weights. 

75 The theory and implementation is inspired from `this paper <http://dx.doi.org/10.1137/S0036144502417715>`_. 

76 

77 Parameters 

78 ---------- 

79 points : list, tuple or np.1darray 

80 The given interpolation points, no specific scaling, but must be 

81 ordered in increasing order. 

82 

83 Attributes 

84 ---------- 

85 points : np.1darray 

86 The interpolating points 

87 weights : np.1darray 

88 The associated barycentric weights 

89 """ 

90 

91 def __init__(self, points, fValues=None): 

92 points = np.asarray(points).ravel() 

93 

94 diffs = points[:, None] - points[None, :] 

95 diffs[np.diag_indices_from(diffs)] = 1 

96 

97 def analytic(diffs): 

98 # Fast implementation (unstable for large number of points) 

99 invProd = np.prod(diffs, axis=1) 

100 invProd **= -1 

101 return invProd 

102 

103 with np.errstate(divide='raise', over='ignore'): 

104 try: 

105 weights = analytic(diffs) 

106 except FloatingPointError: 

107 raise ValueError('Lagrange formula unstable for that much nodes') 

108 

109 # Store attributes 

110 self.points = points 

111 self.weights = weights 

112 

113 # Store function values if provided 

114 if fValues is not None: 

115 fValues = np.asarray(fValues) 

116 if fValues.shape != points.shape: 

117 raise ValueError(f'fValues {fValues.shape} has not the correct shape: {points.shape}') 

118 self.fValues = fValues 

119 

120 def __call__(self, t): 

121 assert self.fValues is not None, "cannot evaluate polynomial without fValues" 

122 t = np.asarray(t) 

123 values = self.getInterpolationMatrix(t.ravel()).dot(self.fValues) 

124 values.shape = t.shape 

125 return values 

126 

127 @property 

128 def n(self): 

129 return self.points.size 

130 

131 def getInterpolationMatrix(self, times): 

132 r""" 

133 Compute the interpolation matrix for a given set of discrete "time" 

134 points. 

135 

136 For instance, if we note :math:`\vec{f}` the vector containing the 

137 :math:`f_j=f(t_j)` values, and :math:`(\tau_m)_{0\leq m<M}` 

138 the "time" points where to interpolate. 

139 Then :math:`I[\vec{f}]`, the vector containing the interpolated 

140 :math:`f(\tau_m)` values, can be obtained using : 

141 

142 .. math:: 

143 I[\vec{f}] = P_{Inter} \vec{f}, 

144 

145 where :math:`P_{Inter}` is the interpolation matrix returned by this 

146 method. 

147 

148 Parameters 

149 ---------- 

150 times : list-like or np.1darray 

151 The discrete "time" points where to interpolate the function. 

152 

153 Returns 

154 ------- 

155 PInter : np.2darray(M, n) 

156 The interpolation matrix, with :math:`M` rows (size of the **times** 

157 parameter) and :math:`n` columns. 

158 

159 """ 

160 # Compute difference between times and Lagrange points 

161 times = np.asarray(times) 

162 with np.errstate(divide='ignore'): 

163 iDiff = 1 / (times[:, None] - self.points[None, :]) 

164 

165 # Find evaluated positions that coincide with one Lagrange point 

166 concom = (iDiff == np.inf) | (iDiff == -np.inf) 

167 i, j = np.where(concom) 

168 

169 # Replace iDiff by on on those lines to get a simple copy of the value 

170 iDiff[i, :] = concom[i, :] 

171 

172 # Compute interpolation matrix using weights 

173 PInter = iDiff * self.weights 

174 PInter /= PInter.sum(axis=-1)[:, None] 

175 

176 return PInter 

177 

178 def getIntegrationMatrix(self, intervals, numQuad='LEGENDRE_NUMPY'): 

179 r""" 

180 Compute the integration matrix for a given set of intervals. 

181 

182 For instance, if we note :math:`\vec{f}` the vector containing the 

183 :math:`f_j=f(t_j)` values, and 

184 :math:`(\tau_{m,left}, \tau_{m,right})_{0\leq m<M}` the different 

185 interval where the function should be integrated using the barycentric 

186 interpolant polynomial. 

187 Then :math:`\Delta[\vec{f}]`, the vector containing the approximations 

188 of 

189 

190 .. math:: 

191 \int_{\tau_{m,left}}^{\tau_{m,right}} f(t)dt, 

192 

193 can be obtained using : 

194 

195 .. math:: 

196 \Delta[\vec{f}] = P_{Integ} \vec{f}, 

197 

198 where :math:`P_{Integ}` is the interpolation matrix returned by this 

199 method. 

200 

201 Parameters 

202 ---------- 

203 intervals : list of pairs 

204 A list of all integration intervals. 

205 numQuad : str, optional 

206 Quadrature rule used to integrate the interpolant barycentric 

207 polynomial. Can be : 

208 

209 - 'LEGENDRE_NUMPY' : Gauss-Legendre rule from Numpy 

210 - 'LEGENDRE_SCIPY' : Gauss-Legendre rule from Scipy 

211 - 'FEJER' : internaly implemented Fejer-I rule 

212 

213 The default is 'LEGENDRE_NUMPY'. 

214 

215 Returns 

216 ------- 

217 PInter : np.2darray(M, n) 

218 The integration matrix, with :math:`M` rows (number of intervals) 

219 and :math:`n` columns. 

220 """ 

221 if numQuad == 'LEGENDRE_NUMPY': 

222 # Legendre gauss rule, integrate exactly polynomials of deg. (2n-1) 

223 iNodes, iWeights = np.polynomial.legendre.leggauss((self.n + 1) // 2) 

224 elif numQuad == 'LEGENDRE_SCIPY': 

225 # Using Legendre scipy implementation 

226 iNodes, iWeights = roots_legendre((self.n + 1) // 2) 

227 elif numQuad == 'FEJER': 

228 # Fejer-I rule, integrate exactly polynomial of deg. n-1 

229 iNodes, iWeights = computeFejerRule(self.n - ((self.n + 1) % 2)) 

230 else: 

231 raise NotImplementedError(f'numQuad={numQuad}') 

232 

233 # Compute quadrature nodes for each interval 

234 intervals = np.array(intervals) 

235 aj, bj = intervals[:, 0][:, None], intervals[:, 1][:, None] 

236 tau, omega = iNodes[None, :], iWeights[None, :] 

237 tEval = (bj - aj) / 2 * tau + (bj + aj) / 2 

238 

239 # Compute the integrand function on nodes 

240 integrand = self.getInterpolationMatrix(tEval.ravel()).T.reshape((-1,) + tEval.shape) 

241 

242 # Apply quadrature rule to integrate 

243 integrand *= omega 

244 integrand *= (bj - aj) / 2 

245 PInter = integrand.sum(axis=-1).T 

246 

247 return PInter