Coverage for pygeodesy/auxilats/auxDST.py: 97%

106 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-23 12:10 -0400

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Discrete Sine Transforms (AuxDST) in Python, transcoded from I{Karney}'s C++ class 

5U{DST<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>} 

6in I{GeographicLib version 2.2+}. 

7 

8Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed 

9under the MIT/X11 License. For more information, see the U{GeographicLib 

10<https://GeographicLib.SourceForge.io>} documentation. 

11 

12@note: Class L{AuxDST} requires package U{numpy<https://PyPI.org/project/numpy>} to be 

13 installed, version 1.16 or newer and needed for I{exact} area calculations. 

14''' 

15# make sure int/int division yields float quotient, see .basics 

16from __future__ import division as _; del _ # PYCHOK semicolon 

17 

18from pygeodesy.basics import isodd, map2, neg, _reverange, _xnumpy 

19from pygeodesy.constants import PI_2, PI_4, isfinite, \ 

20 _0_0, _0_5, _1_0, _naninf 

21# from pygeodesy.fsums import Fsum # from .karney 

22from pygeodesy.karney import _2cos2x, _ALL_DOCS, Fsum, property_RO 

23# from pygeodesy.lazily import _ALL_DOCS # from .karney 

24# from pygeodesy.props import property_RO # from .karney 

25 

26__all__ = () 

27__version__ = '23.08.20' 

28 

29 

30class AuxDST(object): 

31 '''Discrete Sine Transforms (DST) for I{Auxiliary Latitudes}. 

32 

33 @see: I{Karney}'s C++ class U{DST 

34 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>}. 

35 ''' 

36 _N = 0 

37 

38 def __init__(self, N): 

39 '''New L{AuxDST} instance. 

40 

41 @arg N: Size, number of points (C{int}). 

42 ''' 

43 if N > 0: 

44 self._N = int(N) 

45 # kissfft(N, False) # size, inverse 

46 

47 @staticmethod 

48 def evaluate(sinx, cosx, F, *N): 

49 '''Compute the Fourier sum given the sine and cosine of the angle, 

50 using I{Clenshaw} summation C{sum(B{F}[i] * sin((2*i+1) * x))} 

51 for C{i in range(min(len(B{F}), *B{N}))}. 

52 

53 @arg sinx: The sin(I{sigma}) (C{float}). 

54 @arg cosx: The cos(I{sigma}) (C{float}). 

55 @arg F: The Fourier coefficients (C{float}[]). 

56 @arg N: Optional, (smaller) number of terms to evaluate (C{int}). 

57 

58 @return: Precison I{Clenshaw} sum (C{float}). 

59 

60 @see: Methods C{AuxDST.integral} and C{AuxDST.integral2}. 

61 ''' 

62 a = -_2cos2x(cosx, sinx) 

63 if isfinite(a): 

64 Y0, Y1 = Fsum(), Fsum() 

65 n = _len_N(F, *N) 

66 Fn = list(F[:n]) 

67 _F = Fn.pop 

68 if isodd(n): 

69 Y0 -= _F() 

70 while Fn: # Y0, Y1 negated 

71 Y1 -= Y0 * a + _F() 

72 Y0 -= Y1 * a + _F() 

73 r = -float(_Ys(Y0, -Y1, sinx)) 

74 else: 

75 r = _naninf(-a) 

76 return r 

77 

78 @property_RO 

79 def _fft_numpy(self): 

80 '''(INTERNAL) Get the C{numpy.fft} module, I{once}. 

81 ''' 

82 AuxDST._fft_numpy = fft = _xnumpy(AuxDST, 1, 16).fft # overwrite property_RO 

83 return fft 

84 

85 def _fft_real(self, data): 

86 '''(INTERNAL) NumPy's I{kissfft}-like C{transform_real} function, 

87 taking C{float}[:N] B{C{data}} and returning C{complex}[:N*2]. 

88 ''' 

89 # <https://GitHub.com/mborgerding/kissfft/blob/master/test/testkiss.py> 

90 return self._fft_numpy.rfftn(data) 

91 

92 def _ffts(self, data, cIV): 

93 '''(INTERNAL) Compute the DST-III or DST-IV FFTransforms. 

94 

95 @arg data: Elements DST-III[0:N+1] or DST-IV[0:N] (C{float}[]) 

96 with DST_III[0] = 0. 

97 @arg cIV: If C{True} DST-IV, otherwise DST-III. 

98 

99 @return: FFTransforms (C{float}[0:N]). 

100 ''' 

101 t, N = (), self.N 

102 if N > 0: 

103 N2 = N * 2 

104 d = list(data) 

105 # assert len(d) == N + (0 if cIV else 1) 

106 

107 if cIV: # DST-IV 

108 from cmath import exp as _cexp 

109 

110 def _cF(c, j, d=-PI_4 / N): 

111 return c * _cexp(complex(0, d * j)) 

112 

113 i = 0 

114 else: # DST-III 

115 i = 1 

116 # assert d[0] == _0_0 

117 

118 def _cF(c, unused): # PYCHOK redef 

119 return c 

120 

121 d += list(reversed(d[i:N])) # i == len(d) - N 

122 d += list(map(neg, d[:N2])) 

123 c = self._fft_real(d) # complex[0:N*2] 

124 n2 = float(-N2) 

125 t = tuple(_cF(c[j], j).imag / n2 for j in range(1, N2, 2)) 

126 return t 

127 

128 def _ffts2(self, data, F): 

129 '''(INTERNAL) Doubled FFTransforms. 

130 

131 @arg data: Grid centers (C{float}[N]). 

132 @arg F: The transforms (C{float}[N]) 

133 

134 @return: Doubled FFTransforms (C{float}[N*2]). 

135 ''' 

136 def _dmF_2(d, F): 

137 return (d - F) * _0_5 

138 

139 def _dpF_2(d, F): 

140 return (d + F) * _0_5 

141 

142 # Copy DST-IV order N transform to d[0:N] 

143 d = self._ffts(data, True) 

144 N = self._N 

145 # assert len(d) >= N and len(F) >= N 

146 # (DST-IV order N - DST-III order N) / 2 

147 m = map2(_dmF_2, d[:N], F[:N]) 

148 # (DST-IV order N + DST-III order N) / 2 

149 p = map2(_dpF_2, d[:N], F[:N]) 

150 return p + tuple(reversed(m)) 

151 

152 @staticmethod 

153 def integral(sinx, cosx, F, *N): 

154 '''Compute the integral of Fourier sum given the sine and cosine 

155 of the angle using I{Clenshaw} summation C{-sum(B{F}[i] / (2*i+1) * 

156 cos((2*i+1) * x))} for C{i in range(min(len(B{F}), *B{N}))}. 

157 

158 @arg sinx: The sin(I{sigma}) (C{float}). 

159 @arg cosx: The cos(I{sigma}) (C{float}). 

160 @arg F: The Fourier coefficients (C{float}[]). 

161 @arg N: Optional, (smaller) number of terms to evaluate (C{int}). 

162 

163 @return: Precison I{Clenshaw} intergral (C{float}). 

164 

165 @see: Methods C{AuxDST.evaluate} and C{AuxDST.integral2}. 

166 ''' 

167 a = _2cos2x(cosx - sinx, cosx + sinx) 

168 if isfinite(a): 

169 Y0, Y1 = Fsum(), Fsum() 

170 for r in _reverscaled(F, *N): 

171 Y1 -= Y0 * a + r 

172 Y1, Y0 = Y0, -Y1 

173 r = float(_Ys(Y1, Y0, cosx)) 

174 else: 

175 r = _naninf(a) 

176 return r 

177 

178 @staticmethod 

179 def integral2(sin1, cos1, sin2, cos2, F, *N): # PYCHOK no cover 

180 '''Compute the integral of Fourier sum given the sine and cosine 

181 of the angles at the end points using I{Clenshaw} summation 

182 C{integral(siny, cosy, F) - integral(sinx, cosx, F)}. 

183 

184 @arg sin1: The sin(I{sigma1}) (C{float}). 

185 @arg cos1: The cos(I{sigma1}) (C{float}). 

186 @arg sin2: The sin(I{sigma2}) (C{float}). 

187 @arg cos2: The cos(I{sigma2}) (C{float}). 

188 @arg F: The Fourier coefficients (C{float}[]). 

189 @arg N: Optional, (smaller) number of terms to evaluate (C{int}). 

190 

191 @return: Precison I{Clenshaw} intergral (C{float}). 

192 

193 @see: Methods C{AuxDST.evaluate} and C{AuxDST.integral}. 

194 ''' 

195 # 2 * cos(y - x)*cos(y + x) -> 2 * cos(2 * x) 

196 a = _2cos2x(cos2 * cos1, sin2 * sin1) 

197 # -2 * sin(y - x)*sin(y + x) -> 0 

198 b = -_2cos2x(sin2 * cos1, cos2 * sin1) 

199 if isfinite(a) and isfinite(b): 

200 Y0, Y1 = Fsum(), Fsum() 

201 Z0, Z1 = Fsum(), Fsum() 

202 for r in _reverscaled(F, *N): 

203 Y1 -= Y0 * a + Z0 * b + r 

204 Z1 -= Y0 * b + Z0 * a 

205 Y1, Y0 = Y0, -Y1 

206 Z1, Z0 = Z0, -Z1 

207 r = float(_Ys(Y1, Y0, cos2 - cos1) + 

208 _Ys(Z1, Z0, cos2 + cos1)) 

209 else: 

210 r = _naninf(a, b) 

211 return r 

212 

213 @property_RO 

214 def N(self): 

215 '''Get this DST's size, number of points (C{int}). 

216 ''' 

217 return self._N 

218 

219 def refine(self, f, F): 

220 '''Double the number of sampled points on a Fourier series. 

221 

222 @arg f: Single-argument function (C{callable(sigma)} with 

223 C{sigma = PI_4 * j / N for j in range(1, N*2, 2)}. 

224 @arg F: The initial Fourier series coefficients (C{float}[:N]). 

225 

226 @return: Fourier series coefficients (C{float}[:N*2]). 

227 ''' 

228 def _data(_f, N): # [:N] 

229 if N > 0: 

230 d = PI_4 / N 

231 for j in range(1, N*2, 2): 

232 yield _f(d * j) 

233 

234 return self._ffts2(_data(f, self.N), F) 

235 

236 def reset(self, N): 

237 '''Reset this DST. 

238 

239 @arg N: Size, number of points (C{int}). 

240 

241 @return: The new size (C{int}, non-negative). 

242 ''' 

243 self._N = N = max(0, N) 

244 # kissfft.assign(N*2, False) # "reset" size, inverse 

245 return N 

246 

247 def transform(self, f): 

248 '''Compute C{N + 1} terms in the Fourier series. 

249 

250 @arg f: Single-argument function (C{callable(sigma)} with 

251 C{sigma = PI_2 * i / N for i in range(1, N + 1)}. 

252 

253 @return: Fourier series coefficients (C{float}[:N + 1]). 

254 ''' 

255 def _data(_f, N): # [:N + 1] 

256 yield _0_0 # data[0] = 0 

257 if N > 0: 

258 d = PI_2 / N 

259 for i in range(1, N + 1): 

260 yield _f(d * i) 

261 

262 return self._ffts(_data(f, self.N), False) 

263 

264 

265def _len_N(F, *N): 

266 # Adjusted C{len(B{F})}. 

267 return min(len(F), *N) if N else len(F) 

268 

269 

270def _reverscaled(F, *N): 

271 # Yield F[:N], reversed and scaled 

272 for n in _reverange(_len_N(F, *N)): 

273 yield F[n] / (n * 2 + _1_0) 

274 

275 

276def _Ys(X, Y, s): 

277 # Return M{(X - Y) * s}, overwriting X 

278 X -= Y 

279 X *= s 

280 return X 

281 

282 

283__all__ += _ALL_DOCS(AuxDST) 

284 

285# **) MIT License 

286# 

287# Copyright (C) 2023-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

288# 

289# Permission is hereby granted, free of charge, to any person obtaining a 

290# copy of this software and associated documentation files (the "Software"), 

291# to deal in the Software without restriction, including without limitation 

292# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

293# and/or sell copies of the Software, and to permit persons to whom the 

294# Software is furnished to do so, subject to the following conditions: 

295# 

296# The above copyright notice and this permission notice shall be included 

297# in all copies or substantial portions of the Software. 

298# 

299# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

300# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

301# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

302# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

303# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

304# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

305# OTHER DEALINGS IN THE SOFTWARE.