1import numpy as np
2from numpy.testing import assert_allclose
3import pytest
4
5from scipy.fft._fftlog import fht, ifht, fhtoffset
6from scipy.special import poch
7
8
9def test_fht_agrees_with_fftlog():
10    # check that fht numerically agrees with the output from Fortran FFTLog,
11    # the results were generated with the provided `fftlogtest` program,
12    # after fixing how the k array is generated (divide range by n-1, not n)
13
14    # test function, analytical Hankel transform is of the same form
15    def f(r, mu):
16        return r**(mu+1)*np.exp(-r**2/2)
17
18    r = np.logspace(-4, 4, 16)
19
20    dln = np.log(r[1]/r[0])
21    mu = 0.3
22    offset = 0.0
23    bias = 0.0
24
25    a = f(r, mu)
26
27    # test 1: compute as given
28    ours = fht(a, dln, mu, offset=offset, bias=bias)
29    theirs = [ -0.1159922613593045E-02,  0.1625822618458832E-02,
30               -0.1949518286432330E-02,  0.3789220182554077E-02,
31                0.5093959119952945E-03,  0.2785387803618774E-01,
32                0.9944952700848897E-01,  0.4599202164586588    ,
33                0.3157462160881342    , -0.8201236844404755E-03,
34               -0.7834031308271878E-03,  0.3931444945110708E-03,
35               -0.2697710625194777E-03,  0.3568398050238820E-03,
36               -0.5554454827797206E-03,  0.8286331026468585E-03 ]
37    assert_allclose(ours, theirs)
38
39    # test 2: change to optimal offset
40    offset = fhtoffset(dln, mu, bias=bias)
41    ours = fht(a, dln, mu, offset=offset, bias=bias)
42    theirs = [  0.4353768523152057E-04, -0.9197045663594285E-05,
43                0.3150140927838524E-03,  0.9149121960963704E-03,
44                0.5808089753959363E-02,  0.2548065256377240E-01,
45                0.1339477692089897    ,  0.4821530509479356    ,
46                0.2659899781579785    , -0.1116475278448113E-01,
47                0.1791441617592385E-02, -0.4181810476548056E-03,
48                0.1314963536765343E-03, -0.5422057743066297E-04,
49                0.3208681804170443E-04, -0.2696849476008234E-04 ]
50    assert_allclose(ours, theirs)
51
52    # test 3: positive bias
53    bias = 0.8
54    offset = fhtoffset(dln, mu, bias=bias)
55    ours = fht(a, dln, mu, offset=offset, bias=bias)
56    theirs = [ -7.343667355831685     ,  0.1710271207817100    ,
57                0.1065374386206564    , -0.5121739602708132E-01,
58                0.2636649319269470E-01,  0.1697209218849693E-01,
59                0.1250215614723183    ,  0.4739583261486729    ,
60                0.2841149874912028    , -0.8312764741645729E-02,
61                0.1024233505508988E-02, -0.1644902767389120E-03,
62                0.3305775476926270E-04, -0.7786993194882709E-05,
63                0.1962258449520547E-05, -0.8977895734909250E-06 ]
64    assert_allclose(ours, theirs)
65
66    # test 4: negative bias
67    bias = -0.8
68    offset = fhtoffset(dln, mu, bias=bias)
69    ours = fht(a, dln, mu, offset=offset, bias=bias)
70    theirs = [  0.8985777068568745E-05,  0.4074898209936099E-04,
71                0.2123969254700955E-03,  0.1009558244834628E-02,
72                0.5131386375222176E-02,  0.2461678673516286E-01,
73                0.1235812845384476    ,  0.4719570096404403    ,
74                0.2893487490631317    , -0.1686570611318716E-01,
75                0.2231398155172505E-01, -0.1480742256379873E-01,
76                0.1692387813500801    ,  0.3097490354365797    ,
77                2.759360718240186     , 10.52510750700458       ]
78    assert_allclose(ours, theirs)
79
80
81@pytest.mark.parametrize('optimal', [True, False])
82@pytest.mark.parametrize('offset', [0.0, 1.0, -1.0])
83@pytest.mark.parametrize('bias', [0, 0.1, -0.1])
84@pytest.mark.parametrize('n', [64, 63])
85def test_fht_identity(n, bias, offset, optimal):
86    rng = np.random.RandomState(3491349965)
87
88    a = rng.standard_normal(n)
89    dln = rng.uniform(-1, 1)
90    mu = rng.uniform(-2, 2)
91
92    if optimal:
93        offset = fhtoffset(dln, mu, initial=offset, bias=bias)
94
95    A = fht(a, dln, mu, offset=offset, bias=bias)
96    a_ = ifht(A, dln, mu, offset=offset, bias=bias)
97
98    assert_allclose(a, a_)
99
100
101def test_fht_special_cases():
102    rng = np.random.RandomState(3491349965)
103
104    a = rng.standard_normal(64)
105    dln = rng.uniform(-1, 1)
106
107    # let xp = (mu+1+q)/2, xm = (mu+1-q)/2, M = {0, -1, -2, ...}
108
109    # case 1: xp in M, xm in M => well-defined transform
110    mu, bias = -4.0, 1.0
111    with pytest.warns(None) as record:
112        fht(a, dln, mu, bias=bias)
113        assert not record, 'fht warned about a well-defined transform'
114
115    # case 2: xp not in M, xm in M => well-defined transform
116    mu, bias = -2.5, 0.5
117    with pytest.warns(None) as record:
118        fht(a, dln, mu, bias=bias)
119        assert not record, 'fht warned about a well-defined transform'
120
121    # case 3: xp in M, xm not in M => singular transform
122    mu, bias = -3.5, 0.5
123    with pytest.warns(Warning) as record:
124        fht(a, dln, mu, bias=bias)
125        assert record, 'fht did not warn about a singular transform'
126
127    # case 4: xp not in M, xm in M => singular inverse transform
128    mu, bias = -2.5, 0.5
129    with pytest.warns(Warning) as record:
130        ifht(a, dln, mu, bias=bias)
131        assert record, 'ifht did not warn about a singular transform'
132
133
134@pytest.mark.parametrize('n', [64, 63])
135def test_fht_exact(n):
136    rng = np.random.RandomState(3491349965)
137
138    # for a(r) a power law r^\gamma, the fast Hankel transform produces the
139    # exact continuous Hankel transform if biased with q = \gamma
140
141    mu = rng.uniform(0, 3)
142
143    # convergence of HT: -1-mu < gamma < 1/2
144    gamma = rng.uniform(-1-mu, 1/2)
145
146    r = np.logspace(-2, 2, n)
147    a = r**gamma
148
149    dln = np.log(r[1]/r[0])
150
151    offset = fhtoffset(dln, mu, initial=0.0, bias=gamma)
152
153    A = fht(a, dln, mu, offset=offset, bias=gamma)
154
155    k = np.exp(offset)/r[::-1]
156
157    # analytical result
158    At = (2/k)**gamma * poch((mu+1-gamma)/2, gamma)
159
160    assert_allclose(A, At)
161