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)