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

100 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-06 15:27 -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:Charles@Karney.com>} (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.auxilats.auxily import _2cos2x 

19from pygeodesy.basics import isodd, map2, neg, _xnumpy 

20from pygeodesy.constants import PI_2, PI_4, _0_0, _0_5 

21from pygeodesy.fsums import Fsum, Property_RO, property_RO 

22from pygeodesy.lazily import _ALL_DOCS 

23# from pygeodesy.props import Property_RO, property_RO # from .fsums 

24 

25__all__ = () 

26__version__ = '23.08.06' 

27 

28 

29class AuxDST(object): 

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

31 

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

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

34 ''' 

35 _N = 0 

36 

37 def __init__(self, N): 

38 '''New L{AuxDST} instance. 

39 

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

41 ''' 

42 if N > 0: 

43 self._N = int(N) 

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

45 

46 @staticmethod 

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

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

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

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

51 

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

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

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

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

56 

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

58 

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

60 ''' 

61 n, Y1 = len(F), Fsum() 

62 if N: 

63 n = min(n, *N) 

64 if isodd(n): 

65 n -= 1 

66 Y0 = Fsum(-F[n]) 

67 else: 

68 Y0 = Fsum() 

69 a = -_2cos2x(cosx, sinx) 

70 while n > 0: # Y0, Y1 negated 

71 n -= 1; Y1 -= Y0 * a + F[n] # PYCHOK semicolon 

72 n -= 1; Y0 -= Y1 * a + F[n] # PYCHOK semicolon 

73 return -float((Y0 + Y1) * sinx) 

74 

75 @Property_RO 

76 def _fft_real(self): 

77 '''(INTERNAL) Get NumPy's I{kiss}-like C{transform_real}, 

78 taking C{float}[:N] and returning C{complex}[:N*2]. 

79 ''' 

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

81 return _xnumpy(AuxDST, 1, 16).fft.rfftn 

82 

83 def _ffts(self, data, cIV): 

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

85 

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

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

88 

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

90 ''' 

91 N = self.N 

92 if N > 0: 

93 N2 = N * 2 

94 d = list(data) 

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

96 

97 if cIV: # DST-IV 

98 from cmath import exp as _cexp 

99 

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

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

102 

103 i = 0 

104 else: # DST-III 

105 i = 1 

106 # assert d[0] == _0_0 

107 

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

109 return c 

110 

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

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

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

114 n2 = float(-N2) 

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

116 else: 

117 d = () 

118 return d 

119 

120 def _ffts2(self, data, F): 

121 '''(INTERNAL) Doubled FFTransforms. 

122 

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

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

125 

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

127 ''' 

128 def _dmF(d, F): 

129 return (d - F) * _0_5 

130 

131 def _dpF(d, F): 

132 return (d + F) * _0_5 

133 

134 N = self._N 

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

136 d = self._ffts(data, True) 

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

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

139 M = map2(_dmF, d[:N], F[:N]) 

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

141 P = map2(_dpF, d[:N], F[:N]) 

142 return P + tuple(reversed(M)) 

143 

144 @staticmethod 

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

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

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

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

149 

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

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

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

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

154 

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

156 

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

158 ''' 

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

160 Y0, Y1 = Fsum(), Fsum() 

161 for r in _reveN(F, *N): 

162 Y1 -= Y0 * a + r 

163 Y0, Y1 = -Y1, Y0 

164 return float((Y1 - Y0) * cosx) 

165 

166# @staticmethod 

167# def integral2(sin1, cos1, sin2, cos2, F, *N): 

168# '''Compute the integral of Fourier sum given the sine and cosine 

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

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

171# 

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

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

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

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

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

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

178# 

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

180# 

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

182# ''' 

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

184# a = _2cos2x(cos2 * cos1, sin2 * sin1) 

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

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

187# Y0, Y1 = Fsum(), Fsum() 

188# Z0, Z1 = Fsum(), Fsum() 

189# for r in _reveN(F, *N): 

190# Y1 -= Y0 * a + Z0 * b + r 

191# Z1 -= Y0 * b + Z0 * a 

192# Y0, Y1 = -Y1, Y0 

193# Z0, Z1 = -Z1, Z0 

194# return float((Y1 - Y0) * (cos2 - cos1) + 

195# (Z1 - Z0) * (cos2 + cos1)) 

196 

197 @property_RO 

198 def N(self): 

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

200 ''' 

201 return self._N 

202 

203 def refine(self, f, F): 

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

205 

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

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

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

209 

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

211 ''' 

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

213 if N > 0: 

214 d = PI_4 / N 

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

216 yield _f(d * j) 

217 

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

219 

220 def reset(self, N): 

221 '''Reset this DST. 

222 

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

224 

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

226 ''' 

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

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

229 return N 

230 

231 def transform(self, f): 

232 '''Compute C{N} terms in the Fourier series. 

233 

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

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

236 

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

238 ''' 

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

240 yield _0_0 # data[0] = 0 

241 if N > 0: 

242 d = PI_2 / N 

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

244 yield _f(d * i) 

245 

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

247 

248 

249def _reveN(F, *N): 

250 # Yield F reversed and scaled 

251 n = len(F) 

252 if N: 

253 n = min(n, *N) 

254 if n > 0: 

255 n2 = n * 2 + 1 

256 for r in reversed(F[:n]): 

257 n2 -= 2 

258 yield r / n2 

259 

260 

261__all__ += _ALL_DOCS(AuxDST) 

262 

263# **) MIT License 

264# 

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

266# 

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

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

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

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

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

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

273# 

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

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

276# 

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

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

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

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

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

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

283# OTHER DEALINGS IN THE SOFTWARE.