fft_stress.py 6.79 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 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/>.
#
# Copyright(C) 2019-2020 Max-Planck-Society


Martin Reinecke's avatar
Martin Reinecke committed
17
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
18
import ducc0.fft as fft
Martin Reinecke's avatar
Martin Reinecke committed
19
20


21
22
23
rng = np.random.default_rng(42)


Martin Reinecke's avatar
Martin Reinecke committed
24
def _l2error(a, b, axes):
Martin Reinecke's avatar
Martin Reinecke committed
25
    return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))/np.log2(np.max([2, np.prod(np.take(a.shape, axes))]))
Martin Reinecke's avatar
Martin Reinecke committed
26
27
28


def fftn(a, axes=None, inorm=0, out=None, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
29
30
    return fft.c2c(a, axes=axes, forward=True, inorm=inorm,
                   out=out, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
31
32
33


def ifftn(a, axes=None, inorm=0, out=None, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
34
35
    return fft.c2c(a, axes=axes, forward=False, inorm=inorm,
                   out=out, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
36
37
38


def rfftn(a, axes=None, inorm=0, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
39
40
    return fft.r2c(a, axes=axes, forward=True, inorm=inorm,
                   nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43


def irfftn(a, axes=None, lastsize=0, inorm=0, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
44
45
    return fft.c2r(a, axes=axes, lastsize=lastsize, forward=False,
                   inorm=inorm, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
46
47
48
49
50


nthreads = 0


Martin Reinecke's avatar
Martin Reinecke committed
51
def update_err(err, name, value, shape):
Martin Reinecke's avatar
Martin Reinecke committed
52
53
54
    if name in err and err[name] >= value:
        return err
    err[name] = value
Martin Reinecke's avatar
Martin Reinecke committed
55
    print(shape)
Martin Reinecke's avatar
Martin Reinecke committed
56
57
58
59
60
61
62
    for (nm, v) in err.items():
        print("{}: {}".format(nm, v))
    print()
    return err


def test(err):
63
    ndim = rng.integers(1, 5)
Martin Reinecke's avatar
Martin Reinecke committed
64
    axlen = int((2**20)**(1./ndim))
65
    shape = rng.integers(1, axlen, ndim)
Martin Reinecke's avatar
Martin Reinecke committed
66
    axes = np.arange(ndim)
67
68
    rng.shuffle(axes)
    nax = rng.integers(1, ndim+1)
Martin Reinecke's avatar
Martin Reinecke committed
69
70
    axes = axes[:nax]
    lastsize = shape[axes[-1]]
71
    a = rng.random(shape)-0.5 + 1j*rng.random(shape)-0.5j
Martin Reinecke's avatar
Martin Reinecke committed
72
73
74
    a_32 = a.astype(np.complex64)
    b = ifftn(fftn(a, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
              nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
75
    err = update_err(err, "cmax", _l2error(a, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
76
77
    b = ifftn(fftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
              nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
78
    err = update_err(err, "cmax", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
79
80
    b = fftn(ifftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
             nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
81
    err = update_err(err, "cmax", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
82
83
    b = ifftn(fftn(a.astype(np.complex64), axes=axes, nthreads=nthreads),
              axes=axes, inorm=2, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
84
    err = update_err(err, "cmaxf", _l2error(a.astype(np.complex64), b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
85
86
    b = irfftn(rfftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
               lastsize=lastsize, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
87
    err = update_err(err, "rmax", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
88
89
    b = irfftn(rfftn(a.real.astype(np.float32), axes=axes, nthreads=nthreads),
               axes=axes, inorm=2, lastsize=lastsize, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
90
    err = update_err(err, "rmaxf", _l2error(a.real.astype(np.float32), b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
91
92
    b = fft.separable_hartley(
        fft.separable_hartley(a.real, axes=axes, nthreads=nthreads),
Martin Reinecke's avatar
Martin Reinecke committed
93
        axes=axes, inorm=2, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
94
    err = update_err(err, "hmax", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
95
96
    b = fft.genuine_hartley(
        fft.genuine_hartley(a.real, axes=axes, nthreads=nthreads),
Martin Reinecke's avatar
Martin Reinecke committed
97
        axes=axes, inorm=2, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
98
    err = update_err(err, "hmax", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
99
100
    b = fft.separable_hartley(
            fft.separable_hartley(
Martin Reinecke's avatar
Martin Reinecke committed
101
102
                a.real.astype(np.float32), axes=axes, nthreads=nthreads),
            axes=axes, inorm=2, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
103
    err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
104
105
    b = fft.genuine_hartley(
            fft.genuine_hartley(a.real.astype(np.float32), axes=axes,
Martin Reinecke's avatar
Martin Reinecke committed
106
                                nthreads=nthreads),
Martin Reinecke's avatar
Martin Reinecke committed
107
            axes=axes, inorm=2, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
108
    err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
109
    if all(a.shape[i] > 1 for i in axes):
Martin Reinecke's avatar
Martin Reinecke committed
110
111
        b = fft.dct(
            fft.dct(a.real, axes=axes, nthreads=nthreads, type=1),
Martin Reinecke's avatar
Martin Reinecke committed
112
            axes=axes, type=1, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
113
        err = update_err(err, "c1max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
114
115
        b = fft.dct(
            fft.dct(a_32.real, axes=axes, nthreads=nthreads, type=1),
Martin Reinecke's avatar
Martin Reinecke committed
116
            axes=axes, type=1, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
117
        err = update_err(err, "c1maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
118
119
    b = fft.dct(
        fft.dct(a.real, axes=axes, nthreads=nthreads, type=2),
Martin Reinecke's avatar
Martin Reinecke committed
120
        axes=axes, type=3, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
121
    err = update_err(err, "c23max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
122
123
    b = fft.dct(
        fft.dct(a_32.real, axes=axes, nthreads=nthreads, type=2),
Martin Reinecke's avatar
Martin Reinecke committed
124
        axes=axes, type=3, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
125
    err = update_err(err, "c23maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
126
127
    b = fft.dct(
        fft.dct(a.real, axes=axes, nthreads=nthreads, type=4),
Martin Reinecke's avatar
Martin Reinecke committed
128
        axes=axes, type=4, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
129
    err = update_err(err, "c4max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
130
131
    b = fft.dct(
        fft.dct(a_32.real, axes=axes, nthreads=nthreads, type=4),
Martin Reinecke's avatar
Martin Reinecke committed
132
        axes=axes, type=4, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
133
    err = update_err(err, "c4maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
134
135
    b = fft.dst(
        fft.dst(a.real, axes=axes, nthreads=nthreads, type=1),
Martin Reinecke's avatar
Martin Reinecke committed
136
        axes=axes, type=1, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
137
    err = update_err(err, "s1max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
138
139
    b = fft.dst(
        fft.dst(a_32.real, axes=axes, nthreads=nthreads, type=1),
Martin Reinecke's avatar
Martin Reinecke committed
140
        axes=axes, type=1, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
141
    err = update_err(err, "s1maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
142
143
    b = fft.dst(
        fft.dst(a.real, axes=axes, nthreads=nthreads, type=2),
Martin Reinecke's avatar
Martin Reinecke committed
144
        axes=axes, type=3, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
145
    err = update_err(err, "s23max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
146
147
    b = fft.dst(
        fft.dst(a_32.real, axes=axes, nthreads=nthreads, type=2),
Martin Reinecke's avatar
Martin Reinecke committed
148
        axes=axes, type=3, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
149
    err = update_err(err, "s23maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
150
151
    b = fft.dst(
        fft.dst(a.real, axes=axes, nthreads=nthreads, type=4),
Martin Reinecke's avatar
Martin Reinecke committed
152
        axes=axes, type=4, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
153
    err = update_err(err, "s4max", _l2error(a.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
154
155
    b = fft.dst(
        fft.dst(a_32.real, axes=axes, nthreads=nthreads, type=4),
Martin Reinecke's avatar
Martin Reinecke committed
156
        axes=axes, type=4, nthreads=nthreads, inorm=2)
Martin Reinecke's avatar
Martin Reinecke committed
157
    err = update_err(err, "s4maxf", _l2error(a_32.real, b, axes), shape)
Martin Reinecke's avatar
Martin Reinecke committed
158
159
160
161
162


err = dict()
while True:
    test(err)