stress.py 6.46 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.pypocketfft as pypocketfft
Martin Reinecke's avatar
Martin Reinecke committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31


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):
    return pypocketfft.c2c(a, axes=axes, forward=True, inorm=inorm,
                           out=out, nthreads=nthreads)


def ifftn(a, axes=None, inorm=0, out=None, nthreads=1):
    return pypocketfft.c2c(a, axes=axes, forward=False, inorm=inorm,
                           out=out, nthreads=nthreads)


def rfftn(a, axes=None, inorm=0, nthreads=1):
    return pypocketfft.r2c(a, axes=axes, forward=True, inorm=inorm,
                           nthreads=nthreads)


def irfftn(a, axes=None, lastsize=0, inorm=0, nthreads=1):
    return pypocketfft.c2r(a, axes=axes, lastsize=lastsize, forward=False,
                           inorm=inorm, nthreads=nthreads)


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
74
    b = pypocketfft.separable_hartley(
        pypocketfft.separable_hartley(a.real, axes=axes, nthreads=nthreads),
        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
78
    b = pypocketfft.genuine_hartley(
        pypocketfft.genuine_hartley(a.real, axes=axes, nthreads=nthreads),
        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
82
83
    b = pypocketfft.separable_hartley(
            pypocketfft.separable_hartley(
                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
87
88
    b = pypocketfft.genuine_hartley(
            pypocketfft.genuine_hartley(a.real.astype(np.float32), axes=axes,
                                        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
91
92
93
    if all(a.shape[i] > 1 for i in axes):
        b = pypocketfft.dct(
            pypocketfft.dct(a.real, axes=axes, nthreads=nthreads, type=1),
            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
97
        b = pypocketfft.dct(
            pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=1),
            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
101
    b = pypocketfft.dct(
        pypocketfft.dct(a.real, axes=axes, nthreads=nthreads, type=2),
        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
105
    b = pypocketfft.dct(
        pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=2),
        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
109
    b = pypocketfft.dct(
        pypocketfft.dct(a.real, axes=axes, nthreads=nthreads, type=4),
        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
113
    b = pypocketfft.dct(
        pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=4),
        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
117
    b = pypocketfft.dst(
        pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=1),
        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
121
    b = pypocketfft.dst(
        pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=1),
        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
125
    b = pypocketfft.dst(
        pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=2),
        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
129
    b = pypocketfft.dst(
        pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=2),
        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
133
    b = pypocketfft.dst(
        pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=4),
        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
137
    b = pypocketfft.dst(
        pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=4),
        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)