Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/fft/_fftlog.py: 14%
76 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1'''Fast Hankel transforms using the FFTLog algorithm.
3The implementation closely follows the Fortran code of Hamilton (2000).
5added: 14/11/2020 Nicolas Tessore <n.tessore@ucl.ac.uk>
6'''
8import numpy as np
9from warnings import warn
10from ._basic import rfft, irfft
11from ..special import loggamma, poch
13__all__ = [
14 'fht', 'ifht',
15 'fhtoffset',
16]
19# constants
20LN_2 = np.log(2)
23def fht(a, dln, mu, offset=0.0, bias=0.0):
24 r'''Compute the fast Hankel transform.
26 Computes the discrete Hankel transform of a logarithmically spaced periodic
27 sequence using the FFTLog algorithm [1]_, [2]_.
29 Parameters
30 ----------
31 a : array_like (..., n)
32 Real periodic input array, uniformly logarithmically spaced. For
33 multidimensional input, the transform is performed over the last axis.
34 dln : float
35 Uniform logarithmic spacing of the input array.
36 mu : float
37 Order of the Hankel transform, any positive or negative real number.
38 offset : float, optional
39 Offset of the uniform logarithmic spacing of the output array.
40 bias : float, optional
41 Exponent of power law bias, any positive or negative real number.
43 Returns
44 -------
45 A : array_like (..., n)
46 The transformed output array, which is real, periodic, uniformly
47 logarithmically spaced, and of the same shape as the input array.
49 See Also
50 --------
51 ifht : The inverse of `fht`.
52 fhtoffset : Return an optimal offset for `fht`.
54 Notes
55 -----
56 This function computes a discrete version of the Hankel transform
58 .. math::
60 A(k) = \int_{0}^{\infty} \! a(r) \, J_\mu(kr) \, k \, dr \;,
62 where :math:`J_\mu` is the Bessel function of order :math:`\mu`. The index
63 :math:`\mu` may be any real number, positive or negative.
65 The input array `a` is a periodic sequence of length :math:`n`, uniformly
66 logarithmically spaced with spacing `dln`,
68 .. math::
70 a_j = a(r_j) \;, \quad
71 r_j = r_c \exp[(j-j_c) \, \mathtt{dln}]
73 centred about the point :math:`r_c`. Note that the central index
74 :math:`j_c = (n-1)/2` is half-integral if :math:`n` is even, so that
75 :math:`r_c` falls between two input elements. Similarly, the output
76 array `A` is a periodic sequence of length :math:`n`, also uniformly
77 logarithmically spaced with spacing `dln`
79 .. math::
81 A_j = A(k_j) \;, \quad
82 k_j = k_c \exp[(j-j_c) \, \mathtt{dln}]
84 centred about the point :math:`k_c`.
86 The centre points :math:`r_c` and :math:`k_c` of the periodic intervals may
87 be chosen arbitrarily, but it would be usual to choose the product
88 :math:`k_c r_c = k_j r_{n-1-j} = k_{n-1-j} r_j` to be unity. This can be
89 changed using the `offset` parameter, which controls the logarithmic offset
90 :math:`\log(k_c) = \mathtt{offset} - \log(r_c)` of the output array.
91 Choosing an optimal value for `offset` may reduce ringing of the discrete
92 Hankel transform.
94 If the `bias` parameter is nonzero, this function computes a discrete
95 version of the biased Hankel transform
97 .. math::
99 A(k) = \int_{0}^{\infty} \! a_q(r) \, (kr)^q \, J_\mu(kr) \, k \, dr
101 where :math:`q` is the value of `bias`, and a power law bias
102 :math:`a_q(r) = a(r) \, (kr)^{-q}` is applied to the input sequence.
103 Biasing the transform can help approximate the continuous transform of
104 :math:`a(r)` if there is a value :math:`q` such that :math:`a_q(r)` is
105 close to a periodic sequence, in which case the resulting :math:`A(k)` will
106 be close to the continuous transform.
108 References
109 ----------
110 .. [1] Talman J. D., 1978, J. Comp. Phys., 29, 35
111 .. [2] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191)
113 Examples
114 --------
116 This example is the adapted version of ``fftlogtest.f`` which is provided
117 in [2]_. It evaluates the integral
119 .. math::
121 \int^\infty_0 r^{\mu+1} \exp(-r^2/2) J_\mu(k, r) k dr
122 = k^{\mu+1} \exp(-k^2/2) .
124 >>> import numpy as np
125 >>> from scipy import fft
126 >>> import matplotlib.pyplot as plt
128 Parameters for the transform.
130 >>> mu = 0.0 # Order mu of Bessel function
131 >>> r = np.logspace(-7, 1, 128) # Input evaluation points
132 >>> dln = np.log(r[1]/r[0]) # Step size
133 >>> offset = fft.fhtoffset(dln, initial=-6*np.log(10), mu=mu)
134 >>> k = np.exp(offset)/r[::-1] # Output evaluation points
136 Define the analytical function.
138 >>> def f(x, mu):
139 ... """Analytical function: x^(mu+1) exp(-x^2/2)."""
140 ... return x**(mu + 1)*np.exp(-x**2/2)
142 Evaluate the function at ``r`` and compute the corresponding values at
143 ``k`` using FFTLog.
145 >>> a_r = f(r, mu)
146 >>> fht = fft.fht(a_r, dln, mu=mu, offset=offset)
148 For this example we can actually compute the analytical response (which in
149 this case is the same as the input function) for comparison and compute the
150 relative error.
152 >>> a_k = f(k, mu)
153 >>> rel_err = abs((fht-a_k)/a_k)
155 Plot the result.
157 >>> figargs = {'sharex': True, 'sharey': True, 'constrained_layout': True}
158 >>> fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), **figargs)
159 >>> ax1.set_title(r'$r^{\mu+1}\ \exp(-r^2/2)$')
160 >>> ax1.loglog(r, a_r, 'k', lw=2)
161 >>> ax1.set_xlabel('r')
162 >>> ax2.set_title(r'$k^{\mu+1} \exp(-k^2/2)$')
163 >>> ax2.loglog(k, a_k, 'k', lw=2, label='Analytical')
164 >>> ax2.loglog(k, fht, 'C3--', lw=2, label='FFTLog')
165 >>> ax2.set_xlabel('k')
166 >>> ax2.legend(loc=3, framealpha=1)
167 >>> ax2.set_ylim([1e-10, 1e1])
168 >>> ax2b = ax2.twinx()
169 >>> ax2b.loglog(k, rel_err, 'C0', label='Rel. Error (-)')
170 >>> ax2b.set_ylabel('Rel. Error (-)', color='C0')
171 >>> ax2b.tick_params(axis='y', labelcolor='C0')
172 >>> ax2b.legend(loc=4, framealpha=1)
173 >>> ax2b.set_ylim([1e-9, 1e-3])
174 >>> plt.show()
176 '''
178 # size of transform
179 n = np.shape(a)[-1]
181 # bias input array
182 if bias != 0:
183 # a_q(r) = a(r) (r/r_c)^{-q}
184 j_c = (n-1)/2
185 j = np.arange(n)
186 a = a * np.exp(-bias*(j - j_c)*dln)
188 # compute FHT coefficients
189 u = fhtcoeff(n, dln, mu, offset=offset, bias=bias)
191 # transform
192 A = _fhtq(a, u)
194 # bias output array
195 if bias != 0:
196 # A(k) = A_q(k) (k/k_c)^{-q} (k_c r_c)^{-q}
197 A *= np.exp(-bias*((j - j_c)*dln + offset))
199 return A
202def ifht(A, dln, mu, offset=0.0, bias=0.0):
203 r'''Compute the inverse fast Hankel transform.
205 Computes the discrete inverse Hankel transform of a logarithmically spaced
206 periodic sequence. This is the inverse operation to `fht`.
208 Parameters
209 ----------
210 A : array_like (..., n)
211 Real periodic input array, uniformly logarithmically spaced. For
212 multidimensional input, the transform is performed over the last axis.
213 dln : float
214 Uniform logarithmic spacing of the input array.
215 mu : float
216 Order of the Hankel transform, any positive or negative real number.
217 offset : float, optional
218 Offset of the uniform logarithmic spacing of the output array.
219 bias : float, optional
220 Exponent of power law bias, any positive or negative real number.
222 Returns
223 -------
224 a : array_like (..., n)
225 The transformed output array, which is real, periodic, uniformly
226 logarithmically spaced, and of the same shape as the input array.
228 See Also
229 --------
230 fht : Definition of the fast Hankel transform.
231 fhtoffset : Return an optimal offset for `ifht`.
233 Notes
234 -----
235 This function computes a discrete version of the Hankel transform
237 .. math::
239 a(r) = \int_{0}^{\infty} \! A(k) \, J_\mu(kr) \, r \, dk \;,
241 where :math:`J_\mu` is the Bessel function of order :math:`\mu`. The index
242 :math:`\mu` may be any real number, positive or negative.
244 See `fht` for further details.
246 '''
248 # size of transform
249 n = np.shape(A)[-1]
251 # bias input array
252 if bias != 0:
253 # A_q(k) = A(k) (k/k_c)^{q} (k_c r_c)^{q}
254 j_c = (n-1)/2
255 j = np.arange(n)
256 A = A * np.exp(bias*((j - j_c)*dln + offset))
258 # compute FHT coefficients
259 u = fhtcoeff(n, dln, mu, offset=offset, bias=bias)
261 # transform
262 a = _fhtq(A, u, inverse=True)
264 # bias output array
265 if bias != 0:
266 # a(r) = a_q(r) (r/r_c)^{q}
267 a /= np.exp(-bias*(j - j_c)*dln)
269 return a
272def fhtcoeff(n, dln, mu, offset=0.0, bias=0.0):
273 '''Compute the coefficient array for a fast Hankel transform.
274 '''
276 lnkr, q = offset, bias
278 # Hankel transform coefficients
279 # u_m = (kr)^{-i 2m pi/(n dlnr)} U_mu(q + i 2m pi/(n dlnr))
280 # with U_mu(x) = 2^x Gamma((mu+1+x)/2)/Gamma((mu+1-x)/2)
281 xp = (mu+1+q)/2
282 xm = (mu+1-q)/2
283 y = np.linspace(0, np.pi*(n//2)/(n*dln), n//2+1)
284 u = np.empty(n//2+1, dtype=complex)
285 v = np.empty(n//2+1, dtype=complex)
286 u.imag[:] = y
287 u.real[:] = xm
288 loggamma(u, out=v)
289 u.real[:] = xp
290 loggamma(u, out=u)
291 y *= 2*(LN_2 - lnkr)
292 u.real -= v.real
293 u.real += LN_2*q
294 u.imag += v.imag
295 u.imag += y
296 np.exp(u, out=u)
298 # fix last coefficient to be real
299 u.imag[-1] = 0
301 # deal with special cases
302 if not np.isfinite(u[0]):
303 # write u_0 = 2^q Gamma(xp)/Gamma(xm) = 2^q poch(xm, xp-xm)
304 # poch() handles special cases for negative integers correctly
305 u[0] = 2**q * poch(xm, xp-xm)
306 # the coefficient may be inf or 0, meaning the transform or the
307 # inverse transform, respectively, is singular
309 return u
312def fhtoffset(dln, mu, initial=0.0, bias=0.0):
313 '''Return optimal offset for a fast Hankel transform.
315 Returns an offset close to `initial` that fulfils the low-ringing
316 condition of [1]_ for the fast Hankel transform `fht` with logarithmic
317 spacing `dln`, order `mu` and bias `bias`.
319 Parameters
320 ----------
321 dln : float
322 Uniform logarithmic spacing of the transform.
323 mu : float
324 Order of the Hankel transform, any positive or negative real number.
325 initial : float, optional
326 Initial value for the offset. Returns the closest value that fulfils
327 the low-ringing condition.
328 bias : float, optional
329 Exponent of power law bias, any positive or negative real number.
331 Returns
332 -------
333 offset : float
334 Optimal offset of the uniform logarithmic spacing of the transform that
335 fulfils a low-ringing condition.
337 See Also
338 --------
339 fht : Definition of the fast Hankel transform.
341 References
342 ----------
343 .. [1] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191)
345 '''
347 lnkr, q = initial, bias
349 xp = (mu+1+q)/2
350 xm = (mu+1-q)/2
351 y = np.pi/(2*dln)
352 zp = loggamma(xp + 1j*y)
353 zm = loggamma(xm + 1j*y)
354 arg = (LN_2 - lnkr)/dln + (zp.imag + zm.imag)/np.pi
355 return lnkr + (arg - np.round(arg))*dln
358def _fhtq(a, u, inverse=False):
359 '''Compute the biased fast Hankel transform.
361 This is the basic FFTLog routine.
362 '''
364 # size of transform
365 n = np.shape(a)[-1]
367 # check for singular transform or singular inverse transform
368 if np.isinf(u[0]) and not inverse:
369 warn('singular transform; consider changing the bias')
370 # fix coefficient to obtain (potentially correct) transform anyway
371 u = u.copy()
372 u[0] = 0
373 elif u[0] == 0 and inverse:
374 warn('singular inverse transform; consider changing the bias')
375 # fix coefficient to obtain (potentially correct) inverse anyway
376 u = u.copy()
377 u[0] = np.inf
379 # biased fast Hankel transform via real FFT
380 A = rfft(a, axis=-1)
381 if not inverse:
382 # forward transform
383 A *= u
384 else:
385 # backward transform
386 A /= u.conj()
387 A = irfft(A, n, axis=-1)
388 A = A[..., ::-1]
390 return A