Coverage for pygeodesy/auxilats/auxDST.py: 95%
100 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -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: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.
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.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
25__all__ = ()
26__version__ = '23.08.06'
29class AuxDST(object):
30 '''Discrete Sine Transforms (DST) for I{Auxiliary Latitudes}.
32 @see: I{Karney}'s C++ class U{DST
33 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>}.
34 '''
35 _N = 0
37 def __init__(self, N):
38 '''New L{AuxDST} instance.
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
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}))}.
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}).
57 @return: Precison I{Clenshaw} sum (C{float}).
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)
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
83 def _ffts(self, data, cIV):
84 '''(INTERNAL) Compute the DST-III or DST-IV FFTransforms.
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.
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)
97 if cIV: # DST-IV
98 from cmath import exp as _cexp
100 def _cF(c, j, d=-PI_4 / N):
101 return c * _cexp(complex(0, d * j))
103 i = 0
104 else: # DST-III
105 i = 1
106 # assert d[0] == _0_0
108 def _cF(c, unused): # PYCHOK redef
109 return c
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
120 def _ffts2(self, data, F):
121 '''(INTERNAL) Doubled FFTransforms.
123 @arg data: Grid centers (C{float}[N]).
124 @arg F: The transforms (C{float}[N])
126 @return: Doubled FFTransforms (C{float}[N*2]).
127 '''
128 def _dmF(d, F):
129 return (d - F) * _0_5
131 def _dpF(d, F):
132 return (d + F) * _0_5
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))
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}))}.
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}).
155 @return: Precison I{Clenshaw} intergral (C{float}).
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)
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))
197 @property_RO
198 def N(self):
199 '''Get this DST's size, number of points (C{int}).
200 '''
201 return self._N
203 def refine(self, f, F):
204 '''Double the number of sampled points on a Fourier series.
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]).
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)
218 return self._ffts2(_data(f, self.N), F)
220 def reset(self, N):
221 '''Reset this DST.
223 @arg N: Size, number of points (C{int}).
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
231 def transform(self, f):
232 '''Compute C{N} terms in the Fourier series.
234 @arg f: Single-argument function (C{callable(sigma)} with
235 C{sigma = i * PI_2 / N for i in range(1, N + 1)}.
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)
246 return self._ffts(_data(f, self.N), False)
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
261__all__ += _ALL_DOCS(AuxDST)
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.