Coverage for pygeodesy/auxilats/auxDST.py: 97%
106 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-28 15:52 -0400
2# -*- coding: utf-8 -*-
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+}.
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.
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
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
26__all__ = ()
27__version__ = '23.08.20'
30class AuxDST(object):
31 '''Discrete Sine Transforms (DST) for I{Auxiliary Latitudes}.
33 @see: I{Karney}'s C++ class U{DST
34 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>}.
35 '''
36 _N = 0
38 def __init__(self, N):
39 '''New L{AuxDST} instance.
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
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}))}.
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}).
58 @return: Precison I{Clenshaw} sum (C{float}).
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
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
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)
92 def _ffts(self, data, cIV):
93 '''(INTERNAL) Compute the DST-III or DST-IV FFTransforms.
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.
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)
107 if cIV: # DST-IV
108 from cmath import exp as _cexp
110 def _cF(c, j, d=-PI_4 / N):
111 return c * _cexp(complex(0, d * j))
113 i = 0
114 else: # DST-III
115 i = 1
116 # assert d[0] == _0_0
118 def _cF(c, unused): # PYCHOK redef
119 return c
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
128 def _ffts2(self, data, F):
129 '''(INTERNAL) Doubled FFTransforms.
131 @arg data: Grid centers (C{float}[N]).
132 @arg F: The transforms (C{float}[N])
134 @return: Doubled FFTransforms (C{float}[N*2]).
135 '''
136 def _dmF_2(d, F):
137 return (d - F) * _0_5
139 def _dpF_2(d, F):
140 return (d + F) * _0_5
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))
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}))}.
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}).
163 @return: Precison I{Clenshaw} intergral (C{float}).
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
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)}.
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}).
191 @return: Precison I{Clenshaw} intergral (C{float}).
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
213 @property_RO
214 def N(self):
215 '''Get this DST's size, number of points (C{int}).
216 '''
217 return self._N
219 def refine(self, f, F):
220 '''Double the number of sampled points on a Fourier series.
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]).
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)
234 return self._ffts2(_data(f, self.N), F)
236 def reset(self, N):
237 '''Reset this DST.
239 @arg N: Size, number of points (C{int}).
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
247 def transform(self, f):
248 '''Compute C{N + 1} terms in the Fourier series.
250 @arg f: Single-argument function (C{callable(sigma)} with
251 C{sigma = PI_2 * i / N for i in range(1, N + 1)}.
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)
262 return self._ffts(_data(f, self.N), False)
265def _len_N(F, *N):
266 # Adjusted C{len(B{F})}.
267 return min(len(F), *N) if N else len(F)
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)
276def _Ys(X, Y, s):
277 # Return M{(X - Y) * s}, overwriting X
278 X -= Y
279 X *= s
280 return X
283__all__ += _ALL_DOCS(AuxDST)
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.