test_fft_operator.py 5.6 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
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
16
17
#
# 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
Martin Reinecke's avatar
Martin Reinecke committed
22
from numpy.testing import assert_allclose
Martin Reinecke's avatar
Martin Reinecke committed
23
import nifty2go as ift
Martin Reinecke's avatar
Martin Reinecke committed
24
25
from itertools import product
from test.common import expand
26

Theo Steininger's avatar
Theo Steininger committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
def _get_rtol(tp):
Theo Steininger's avatar
Theo Steininger committed
29
    if (tp == np.float64) or (tp == np.complex128):
Martin Reinecke's avatar
Martin Reinecke committed
30
31
32
        return 1e-10
    else:
        return 1e-5
Martin Reinecke's avatar
Martin Reinecke committed
33

Theo Steininger's avatar
Theo Steininger committed
34

35
class FFTOperatorTests(unittest.TestCase):
Martin Reinecke's avatar
Martin Reinecke committed
36
    @expand(product([16, ], [0.1, 1, 3.7],
Martin Reinecke's avatar
Martin Reinecke committed
37
38
                    [np.float64, np.float32, np.complex64, np.complex128]))
    def test_fft1D(self, dim1, d, itp):
Theo Steininger's avatar
Theo Steininger committed
39
        tol = _get_rtol(itp)
Martin Reinecke's avatar
Martin Reinecke committed
40
41
42
        a = ift.RGSpace(dim1, distances=d)
        b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
        fft = ift.FFTOperator(domain=a, target=b)
43

44
        np.random.seed(16)
Martin Reinecke's avatar
Martin Reinecke committed
45
46
        inp = ift.Field.from_random(domain=a, random_type='normal',
                                    std=7, mean=3, dtype=itp)
47
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
48
49
        assert_allclose(ift.dobj.to_global_data(inp.val),
                        ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
Martin Reinecke's avatar
Martin Reinecke committed
50

Martin Reinecke's avatar
Martin Reinecke committed
51
    @expand(product([12, 15], [9, 12], [0.1, 1, 3.7],
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
52
                    [0.4, 1, 2.7],
Martin Reinecke's avatar
Martin Reinecke committed
53
54
                    [np.float64, np.float32, np.complex64, np.complex128]))
    def test_fft2D(self, dim1, dim2, d1, d2, itp):
Theo Steininger's avatar
Theo Steininger committed
55
        tol = _get_rtol(itp)
Martin Reinecke's avatar
Martin Reinecke committed
56
57
58
59
        a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
        b = ift.RGSpace([dim1, dim2],
                        distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True)
        fft = ift.FFTOperator(domain=a, target=b)
60

Martin Reinecke's avatar
Martin Reinecke committed
61
62
        inp = ift.Field.from_random(domain=a, random_type='normal',
                                    std=7, mean=3, dtype=itp)
63
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
64
65
        assert_allclose(ift.dobj.to_global_data(inp.val),
                        ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
66

Martin Reinecke's avatar
Martin Reinecke committed
67
    @expand(product([0, 1, 2],
Martin Reinecke's avatar
Martin Reinecke committed
68
69
                    [np.float64, np.float32, np.complex64, np.complex128]))
    def test_composed_fft(self, index, dtype):
70
        tol = _get_rtol(dtype)
Martin Reinecke's avatar
Martin Reinecke committed
71
72
73
        a = [a1, a2, a3] = [ift.RGSpace((32,)), ift.RGSpace((4, 4)),
                            ift.RGSpace((5, 6))]
        fft = ift.FFTOperator(domain=a, space=index)
74

Martin Reinecke's avatar
Martin Reinecke committed
75
76
        inp = ift.Field.from_random(domain=(a1, a2, a3), random_type='normal',
                                    std=7, mean=3, dtype=dtype)
77
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
78
79
        assert_allclose(ift.dobj.to_global_data(inp.val),
                        ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
Martin Reinecke's avatar
Martin Reinecke committed
80

Theo Steininger's avatar
Theo Steininger committed
81
    @expand(product([0, 3, 6, 11, 30],
82
                    [np.float64, np.float32, np.complex64, np.complex128]))
Theo Steininger's avatar
Theo Steininger committed
83
84
    def test_sht(self, lm, tp):
        tol = _get_rtol(tp)
Martin Reinecke's avatar
Martin Reinecke committed
85
86
87
88
89
        a = ift.LMSpace(lmax=lm)
        b = ift.GLSpace(nlat=lm+1)
        fft = ift.FFTOperator(domain=a, target=b)
        inp = ift.Field.from_random(domain=a, random_type='normal',
                                    std=7, mean=3, dtype=tp)
90
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
91
92
        assert_allclose(ift.dobj.to_global_data(inp.val),
                        ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
93

Theo Steininger's avatar
Theo Steininger committed
94
    @expand(product([128, 256],
95
                    [np.float64, np.float32, np.complex64, np.complex128]))
Theo Steininger's avatar
Theo Steininger committed
96
    def test_sht2(self, lm, tp):
Martin Reinecke's avatar
Martin Reinecke committed
97
98
99
100
101
        a = ift.LMSpace(lmax=lm)
        b = ift.HPSpace(nside=lm//2)
        fft = ift.FFTOperator(domain=a, target=b)
        inp = ift.Field.from_random(domain=a, random_type='normal',
                                    std=1, mean=0, dtype=tp)
102
        out = fft.adjoint_times(fft.times(inp))
Martin Reinecke's avatar
Martin Reinecke committed
103
104
        assert_allclose(ift.dobj.to_global_data(inp.val),
                        ift.dobj.to_global_data(out.val), rtol=1e-3, atol=1e-1)
105
106

    @expand(product([128, 256],
107
                    [np.float64, np.float32, np.complex64, np.complex128]))
108
    def test_dotsht(self, lm, tp):
Martin Reinecke's avatar
Martin Reinecke committed
109
        tol = 10 * _get_rtol(tp)
Martin Reinecke's avatar
Martin Reinecke committed
110
111
112
113
114
        a = ift.LMSpace(lmax=lm)
        b = ift.GLSpace(nlat=lm+1)
        fft = ift.FFTOperator(domain=a, target=b)
        inp = ift.Field.from_random(domain=a, random_type='normal',
                                    std=1, mean=0, dtype=tp)
115
        out = fft.times(inp)
Martin Reinecke's avatar
Martin Reinecke committed
116
117
        v1 = np.sqrt(out.vdot(out))
        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
Martin Reinecke's avatar
Martin Reinecke committed
118
        assert_allclose(v1, v2, rtol=tol, atol=tol)
119
120

    @expand(product([128, 256],
121
                    [np.float64, np.float32, np.complex64, np.complex128]))
122
    def test_dotsht2(self, lm, tp):
Martin Reinecke's avatar
Martin Reinecke committed
123
        tol = 10 * _get_rtol(tp)
Martin Reinecke's avatar
Martin Reinecke committed
124
125
126
127
        a = ift.LMSpace(lmax=lm)
        b = ift.HPSpace(nside=lm//2)
        fft = ift.FFTOperator(domain=a, target=b)
        inp = ift.Field.from_random(domain=a, random_type='normal',
Martin Reinecke's avatar
Martin Reinecke committed
128
                                    std=1, mean=0, dtype=tp)
129
        out = fft.times(inp)
Martin Reinecke's avatar
Martin Reinecke committed
130
131
        v1 = np.sqrt(out.vdot(out))
        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
Martin Reinecke's avatar
Martin Reinecke committed
132
        assert_allclose(v1, v2, rtol=tol, atol=tol)