stress.py 6.26 KB
 Martin Reinecke committed Feb 23, 2020 1 2 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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 ``````import numpy as np import pypocketfft 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 def update_err(err, name, value): if name in err and err[name] >= value: return err err[name] = value 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) err = update_err(err, "cmax", _l2error(a, b, axes)) b = ifftn(fftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "cmax", _l2error(a.real, b, axes)) b = fftn(ifftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "cmax", _l2error(a.real, b, axes)) b = ifftn(fftn(a.astype(np.complex64), axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "cmaxf", _l2error(a.astype(np.complex64), b, axes)) b = irfftn(rfftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2, lastsize=lastsize, nthreads=nthreads) err = update_err(err, "rmax", _l2error(a.real, b, axes)) b = irfftn(rfftn(a.real.astype(np.float32), axes=axes, nthreads=nthreads), axes=axes, inorm=2, lastsize=lastsize, nthreads=nthreads) err = update_err(err, "rmaxf", _l2error(a.real.astype(np.float32), b, axes)) b = pypocketfft.separable_hartley( pypocketfft.separable_hartley(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "hmax", _l2error(a.real, b, axes)) b = pypocketfft.genuine_hartley( pypocketfft.genuine_hartley(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "hmax", _l2error(a.real, b, axes)) b = pypocketfft.separable_hartley( pypocketfft.separable_hartley( a.real.astype(np.float32), axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b, axes)) b = pypocketfft.genuine_hartley( pypocketfft.genuine_hartley(a.real.astype(np.float32), axes=axes, nthreads=nthreads), axes=axes, inorm=2, nthreads=nthreads) err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b, axes)) 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) err = update_err(err, "c1max", _l2error(a.real, b, axes)) b = pypocketfft.dct( pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=1), axes=axes, type=1, nthreads=nthreads, inorm=2) err = update_err(err, "c1maxf", _l2error(a_32.real, b, axes)) b = pypocketfft.dct( pypocketfft.dct(a.real, axes=axes, nthreads=nthreads, type=2), axes=axes, type=3, nthreads=nthreads, inorm=2) err = update_err(err, "c23max", _l2error(a.real, b, axes)) b = pypocketfft.dct( pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=2), axes=axes, type=3, nthreads=nthreads, inorm=2) err = update_err(err, "c23maxf", _l2error(a_32.real, b, axes)) b = pypocketfft.dct( pypocketfft.dct(a.real, axes=axes, nthreads=nthreads, type=4), axes=axes, type=4, nthreads=nthreads, inorm=2) err = update_err(err, "c4max", _l2error(a.real, b, axes)) b = pypocketfft.dct( pypocketfft.dct(a_32.real, axes=axes, nthreads=nthreads, type=4), axes=axes, type=4, nthreads=nthreads, inorm=2) err = update_err(err, "c4maxf", _l2error(a_32.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=1), axes=axes, type=1, nthreads=nthreads, inorm=2) err = update_err(err, "s1max", _l2error(a.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=1), axes=axes, type=1, nthreads=nthreads, inorm=2) err = update_err(err, "s1maxf", _l2error(a_32.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=2), axes=axes, type=3, nthreads=nthreads, inorm=2) err = update_err(err, "s23max", _l2error(a.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=2), axes=axes, type=3, nthreads=nthreads, inorm=2) err = update_err(err, "s23maxf", _l2error(a_32.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a.real, axes=axes, nthreads=nthreads, type=4), axes=axes, type=4, nthreads=nthreads, inorm=2) err = update_err(err, "s4max", _l2error(a.real, b, axes)) b = pypocketfft.dst( pypocketfft.dst(a_32.real, axes=axes, nthreads=nthreads, type=4), axes=axes, type=4, nthreads=nthreads, inorm=2) err = update_err(err, "s4maxf", _l2error(a_32.real, b, axes)) err = dict() while True: test(err)``````