test_fft_operator.py 5.69 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
Martin Reinecke's avatar
Martin Reinecke committed
18

Martin Reinecke's avatar
Martin Reinecke committed
19
from __future__ import division
Martin Reinecke's avatar
Martin Reinecke committed
20
21
import unittest
import numpy as np
Theo Steininger's avatar
Theo Steininger committed
22
23
from numpy.testing import assert_equal,\
    assert_allclose
Martin Reinecke's avatar
Martin Reinecke committed
24
from nifty2go import Field,\
Theo Steininger's avatar
Theo Steininger committed
25
26
    RGSpace,\
    LMSpace,\
27
28
    HPSpace,\
    GLSpace,\
29
    FFTOperator
Martin Reinecke's avatar
Martin Reinecke committed
30
31
from itertools import product
from test.common import expand
Martin Reinecke's avatar
Martin Reinecke committed
32
from nifty2go.dobj import to_ndarray as to_np, from_ndarray as from_np
33

Theo Steininger's avatar
Theo Steininger committed
34

Martin Reinecke's avatar
Martin Reinecke committed
35
def _get_rtol(tp):
Theo Steininger's avatar
Theo Steininger committed
36
    if (tp == np.float64) or (tp == np.complex128):
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
        return 1e-10
    else:
        return 1e-5
Martin Reinecke's avatar
Martin Reinecke committed
40

Theo Steininger's avatar
Theo Steininger committed
41

42
class FFTOperatorTests(unittest.TestCase):
Martin Reinecke's avatar
Martin Reinecke committed
43
    @expand(product([16, ], [0.1, 1, 3.7],
44
45
                    [np.float64, np.float32, np.complex64, np.complex128],
                    ['real', 'complex']))
Martin Reinecke's avatar
Martin Reinecke committed
46
    def test_fft1D(self, dim1, d, itp, base):
Theo Steininger's avatar
Theo Steininger committed
47
        tol = _get_rtol(itp)
48
49
        a = RGSpace(dim1, distances=d)
        b = RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
Martin Reinecke's avatar
Martin Reinecke committed
50
        fft = FFTOperator(domain=a, target=b)
51
52
53
        fft._forward_transformation.harmonic_base = base
        fft._backward_transformation.harmonic_base = base

54
        np.random.seed(16)
Theo Steininger's avatar
Theo Steininger committed
55
        inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
56
                                dtype=itp)
57
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
58
        assert_allclose(to_np(inp.val), to_np(out.val), rtol=tol, atol=tol)
Martin Reinecke's avatar
Martin Reinecke committed
59

Martin Reinecke's avatar
Martin Reinecke committed
60
    @expand(product([12, 15], [9, 12], [0.1, 1, 3.7],
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
61
                    [0.4, 1, 2.7],
62
63
                    [np.float64, np.float32, np.complex64, np.complex128],
                    ['real', 'complex']))
Martin Reinecke's avatar
Martin Reinecke committed
64
    def test_fft2D(self, dim1, dim2, d1, d2,
65
                   itp, base):
Theo Steininger's avatar
Theo Steininger committed
66
        tol = _get_rtol(itp)
67
68
        a = RGSpace([dim1, dim2], distances=[d1, d2])
        b = RGSpace([dim1, dim2],
Martin Reinecke's avatar
Martin Reinecke committed
69
                    distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True)
Martin Reinecke's avatar
Martin Reinecke committed
70
        fft = FFTOperator(domain=a, target=b)
71
72
73
        fft._forward_transformation.harmonic_base = base
        fft._backward_transformation.harmonic_base = base

Theo Steininger's avatar
Theo Steininger committed
74
        inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
75
                                dtype=itp)
76
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
77
        assert_allclose(to_np(inp.val), to_np(out.val), rtol=tol, atol=tol)
78

Martin Reinecke's avatar
Martin Reinecke committed
79
    @expand(product([0, 1, 2],
80
81
                    [np.float64, np.float32, np.complex64, np.complex128],
                    ['real', 'complex']))
Martin Reinecke's avatar
Martin Reinecke committed
82
    def test_composed_fft(self, index, dtype,
83
84
                          base):
        tol = _get_rtol(dtype)
85
        a = [a1, a2, a3] = [RGSpace((32,)), RGSpace((4, 4)), RGSpace((5, 6))]
86
        fft = FFTOperator(domain=a, space=index)
87
88

        inp = Field.from_random(domain=(a1, a2, a3), random_type='normal',
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
89
                                std=7, mean=3, dtype=dtype)
90
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
91
        assert_allclose(to_np(inp.val), to_np(out.val), rtol=tol, atol=tol)
Martin Reinecke's avatar
Martin Reinecke committed
92

Theo Steininger's avatar
Theo Steininger committed
93
    @expand(product([0, 3, 6, 11, 30],
94
                    [np.float64, np.float32, np.complex64, np.complex128]))
Theo Steininger's avatar
Theo Steininger committed
95
96
    def test_sht(self, lm, tp):
        tol = _get_rtol(tp)
Martin Reinecke's avatar
Martin Reinecke committed
97
        a = LMSpace(lmax=lm)
98
        b = GLSpace(nlat=lm+1)
99
        fft = FFTOperator(domain=a, target=b)
Theo Steininger's avatar
Theo Steininger committed
100
        inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3,
Martin Reinecke's avatar
Martin Reinecke committed
101
                                dtype=tp)
102
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
103
        assert_allclose(to_np(inp.val), to_np(out.val), rtol=tol, atol=tol)
104

Theo Steininger's avatar
Theo Steininger committed
105
    @expand(product([128, 256],
106
                    [np.float64, np.float32, np.complex64, np.complex128]))
Theo Steininger's avatar
Theo Steininger committed
107
    def test_sht2(self, lm, tp):
Martin Reinecke's avatar
Martin Reinecke committed
108
        a = LMSpace(lmax=lm)
109
        b = HPSpace(nside=lm//2)
110
        fft = FFTOperator(domain=a, target=b)
Theo Steininger's avatar
Theo Steininger committed
111
        inp = Field.from_random(domain=a, random_type='normal', std=1, mean=0,
Martin Reinecke's avatar
Martin Reinecke committed
112
                                dtype=tp)
113
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
114
        assert_allclose(to_np(inp.val), to_np(out.val), rtol=1e-3, atol=1e-1)
115
116

    @expand(product([128, 256],
117
                    [np.float64, np.float32, np.complex64, np.complex128]))
118
119
120
121
    def test_dotsht(self, lm, tp):
        tol = _get_rtol(tp)
        a = LMSpace(lmax=lm)
        b = GLSpace(nlat=lm+1)
122
        fft = FFTOperator(domain=a, target=b)
123
124
125
        inp = Field.from_random(domain=a, random_type='normal', std=1, mean=0,
                                dtype=tp)
        out = fft.times(inp)
Martin Reinecke's avatar
Martin Reinecke committed
126
127
        v1 = np.sqrt(out.vdot(out))
        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
Martin Reinecke's avatar
Martin Reinecke committed
128
        assert_allclose(v1, v2, rtol=tol, atol=tol)
129
130

    @expand(product([128, 256],
131
                    [np.float64, np.float32, np.complex64, np.complex128]))
132
133
134
135
    def test_dotsht2(self, lm, tp):
        tol = _get_rtol(tp)
        a = LMSpace(lmax=lm)
        b = HPSpace(nside=lm//2)
136
        fft = FFTOperator(domain=a, target=b)
137
138
139
        inp = Field.from_random(domain=a, random_type='normal', std=1, mean=0,
                                dtype=tp)
        out = fft.times(inp)
Martin Reinecke's avatar
Martin Reinecke committed
140
141
        v1 = np.sqrt(out.vdot(out))
        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
Martin Reinecke's avatar
Martin Reinecke committed
142
        assert_allclose(v1, v2, rtol=tol, atol=tol)