stress.py 6.13 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
2
import ducc_0_1.fft as fft
Martin Reinecke's avatar
Martin Reinecke committed
3
4
5
6
7
8
9


def _l2error(a, b, axes):
    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))]))


def fftn(a, axes=None, inorm=0, out=None, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
10
11
    return fft.c2c(a, axes=axes, forward=True, inorm=inorm,
                   out=out, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
12
13
14


def ifftn(a, axes=None, inorm=0, out=None, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
15
16
    return fft.c2c(a, axes=axes, forward=False, inorm=inorm,
                   out=out, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
17
18
19


def rfftn(a, axes=None, inorm=0, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
20
21
    return fft.r2c(a, axes=axes, forward=True, inorm=inorm,
                   nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
22
23
24


def irfftn(a, axes=None, lastsize=0, inorm=0, nthreads=1):
Martin Reinecke's avatar
Martin Reinecke committed
25
26
    return fft.c2r(a, axes=axes, lastsize=lastsize, forward=False,
                   inorm=inorm, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
27
28
29
30
31


nthreads = 0


Martin Reinecke's avatar
Martin Reinecke committed
32
def update_err(err, name, value, shape):
Martin Reinecke's avatar
Martin Reinecke committed
33
34
35
    if name in err and err[name] >= value:
        return err
    err[name] = value
Martin Reinecke's avatar
Martin Reinecke committed
36
    print(shape)
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    for (nm, v) in err.items():
        print("{}: {}".format(nm, v))
    print()
    return err


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


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