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