Commit eb9c915c authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup test programs

parent 87ac8357
import numpy as np
import pypocketfft
import pyfftw
from time import time
import matplotlib.pyplot as plt
import math
nthr = 1
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
def prime_factorize(n):
factors = []
number = math.fabs(n)
while number > 1:
factor = get_next_prime_factor(number)
factors.append(factor)
number /= factor
if n < -1: # If we'd check for < 0, -1 would give us trouble
factors[0] = -factors[0]
return tuple(factors)
def get_next_prime_factor(n):
if n % 2 == 0:
return 2
# Not 'good' [also] checking non-prime numbers I guess?
# But the alternative, creating a list of prime numbers,
# wouldn't it be more demanding? Process of creating it.
for x in range(3, int(math.ceil(math.sqrt(n)) + 1), 2):
if n % x == 0:
return x
return int(n)
def measure_fftw(a, nrepeat):
import pyfftw
f1 = pyfftw.empty_aligned(a.shape, dtype=a.dtype)
f1[()] = a
f2 = pyfftw.empty_aligned(a.shape, dtype=a.dtype)
fftw = pyfftw.FFTW(f1, f2, flags=('FFTW_MEASURE',), threads=nthr)
tmin = 1e38
for i in range(nrepeat):
t0 = time()
fftw()
t1 = time()
tmin = min(tmin, t1-t0)
return tmin
def measure_fftw_np_interface(a, nrepeat):
import pyfftw
pyfftw.interfaces.cache.enable()
tmin = 1e38
for i in range(nrepeat):
t0 = time()
b = pyfftw.interfaces.numpy_fft.fftn(a)
t1 = time()
tmin = min(tmin, t1-t0)
return tmin
def measure_pypocketfft(a, nrepeat):
import pypocketfft as ppf
tmin = 1e38
for i in range(nrepeat):
t0 = time()
b = ppf.c2c(a, forward=True, nthreads=nthr)
t1 = time()
tmin = min(tmin, t1-t0)
return tmin
def measure_scipy_fftpack(a, nrepeat):
import scipy.fftpack
tmin = 1e38
for i in range(nrepeat):
t0 = time()
b = scipy.fftpack.fftn(a)
t1 = time()
tmin = min(tmin, t1-t0)
return tmin
def bench_nd(ndim, nmax, ntry, tp, funcs, nrepeat, ttl="", filename=""):
results = [[] for i in range(len(funcs))]
for n in range(ntry):
print(n, ntry)
shp = np.random.randint(nmax//3, nmax+1, ndim)
a = (np.random.rand(*shp) + 1j*np.random.rand(*shp)).astype(tp)
for func, res in zip(funcs, results):
res.append(func(a, nrepeat))
results = np.array(results)
plt.title("{}: {}D, {}, max_extent={}".format(
ttl, ndim, str(tp), nmax))
plt.xlabel("time ratio")
plt.ylabel("counts")
plt.hist(results[0, :]/results[1, :], bins="auto")
if filename != "":
plt.savefig(filename)
plt.show()
funcs = (measure_pypocketfft, measure_fftw_np_interface)
ttl = "pypocketfft/fftw_numpy_interface"
bench_nd(1, 8192, 100, "c16", funcs, 10, ttl, "1d.png")
bench_nd(2, 2048, 100, "c16", funcs, 2, ttl, "2d.png")
bench_nd(3, 256, 100, "c16", funcs, 2, ttl, "3d.png")
bench_nd(1, 8192, 100, "c8", funcs, 10, ttl, "1d_single.png")
bench_nd(2, 2048, 100, "c8", funcs, 2, ttl, "2d_single.png")
bench_nd(3, 256, 100, "c8", funcs, 2, ttl, "3d_single.png")
...@@ -21,9 +21,9 @@ if sys.platform == 'darwin': ...@@ -21,9 +21,9 @@ if sys.platform == 'darwin':
extra_compile_args += ['--stdlib=libc++', '-mmacosx-version-min=10.9'] extra_compile_args += ['--stdlib=libc++', '-mmacosx-version-min=10.9']
vars = distutils.sysconfig.get_config_vars() vars = distutils.sysconfig.get_config_vars()
vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '') vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '')
python_module_link_args+=['-bundle'] python_module_link_args += ['-bundle']
else: else:
extra_compile_args += ['-DPOCKETFFT_OPENMP', '-fopenmp', '-Wfatal-errors', '-Wfloat-conversion' ,'-Wsign-conversion', '-Wconversion' ,'-W', '-Wall', '-Wstrict-aliasing=2', '-Wwrite-strings', '-Wredundant-decls', '-Woverloaded-virtual', '-Wcast-qual', '-Wcast-align', '-Wpointer-arith'] extra_compile_args += ['-DPOCKETFFT_OPENMP', '-fopenmp', '-Wfatal-errors', '-Wfloat-conversion', '-Wsign-conversion', '-Wconversion' ,'-W', '-Wall', '-Wstrict-aliasing=2', '-Wwrite-strings', '-Wredundant-decls', '-Woverloaded-virtual', '-Wcast-qual', '-Wcast-align', '-Wpointer-arith']
python_module_link_args += ['-march=native', '-Wl,-rpath,$ORIGIN', '-fopenmp'] python_module_link_args += ['-march=native', '-Wl,-rpath,$ORIGIN', '-fopenmp']
# if you don't want debugging info, add "-s" to python_module_link_args # if you don't want debugging info, add "-s" to python_module_link_args
......
import numpy as np import numpy as np
import pypocketfft import pypocketfft
def _l2error(a,b):
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
def fftn(a, axes=None, inorm=0, out=None, nthreads=1): def fftn(a, axes=None, inorm=0, out=None, nthreads=1):
return pypocketfft.c2c(a, axes=axes, forward=True, inorm=inorm, return pypocketfft.c2c(a, axes=axes, forward=True, inorm=inorm,
out=out, nthreads=nthreads) out=out, nthreads=nthreads)
def ifftn(a, axes=None, inorm=0, out=None, nthreads=1): def ifftn(a, axes=None, inorm=0, out=None, nthreads=1):
return pypocketfft.c2c(a, axes=axes, forward=False, inorm=inorm, return pypocketfft.c2c(a, axes=axes, forward=False, inorm=inorm,
out=out, nthreads=nthreads) out=out, nthreads=nthreads)
def rfftn(a, axes=None, inorm=0, nthreads=1): def rfftn(a, axes=None, inorm=0, nthreads=1):
return pypocketfft.r2c(a, axes=axes, forward=True, inorm=inorm, return pypocketfft.r2c(a, axes=axes, forward=True, inorm=inorm,
nthreads=nthreads) nthreads=nthreads)
def irfftn(a, axes=None, lastsize=0, inorm=0, nthreads=1): def irfftn(a, axes=None, lastsize=0, inorm=0, nthreads=1):
return pypocketfft.c2r(a, axes=axes, lastsize=lastsize,forward=False, return pypocketfft.c2r(a, axes=axes, lastsize=lastsize, forward=False,
inorm=inorm, nthreads=nthreads) inorm=inorm, nthreads=nthreads)
nthreads=0
cmaxerr=0. nthreads = 0
fmaxerr=0.
cmaxerrf=0.
fmaxerrf=0. def update_err(err, name, value):
hmaxerr=0. if name not in err:
hmaxerrf=0. err[name] = value
def test(): print("{}: {}".format(name, value))
global cmaxerr, fmaxerr, hmaxerr, cmaxerrf, fmaxerrf, hmaxerrf else:
ndim=np.random.randint(1,5) if value > err[name]:
axlen=int((2**20)**(1./ndim)) err[name] = value
shape = np.random.randint(1,axlen,ndim) print("{}: {}".format(name, value))
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) axes = np.arange(ndim)
np.random.shuffle(axes) np.random.shuffle(axes)
nax = np.random.randint(1,ndim+1) nax = np.random.randint(1, ndim+1)
axes = axes[:nax] axes = axes[:nax]
lastsize = shape[axes[-1]] lastsize = shape[axes[-1]]
a=np.random.rand(*shape)-0.5 + 1j*np.random.rand(*shape)-0.5j a = np.random.rand(*shape)-0.5 + 1j*np.random.rand(*shape)-0.5j
b=ifftn(fftn(a,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) b = ifftn(fftn(a, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
err = _l2error(a,b) nthreads=nthreads)
if err > cmaxerr: err = update_err(err, "cmax", _l2error(a, b))
cmaxerr = err b = ifftn(fftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
print("cmaxerr:", cmaxerr, shape, axes) nthreads=nthreads)
b=ifftn(fftn(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = update_err(err, "cmax", _l2error(a.real, b))
err = _l2error(a.real,b) b = fftn(ifftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
if err > cmaxerr: nthreads=nthreads)
cmaxerr = err err = update_err(err, "cmax", _l2error(a.real, b))
print("cmaxerr:", cmaxerr, shape, axes) b = irfftn(rfftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
b=fftn(ifftn(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) lastsize=lastsize, nthreads=nthreads)
err = _l2error(a.real,b) err = update_err(err, "rmax", _l2error(a.real, b))
if err > cmaxerr: b = ifftn(fftn(a.astype(np.complex64), axes=axes, nthreads=nthreads),
cmaxerr = err axes=axes, inorm=2, nthreads=nthreads)
print("cmaxerr:", cmaxerr, shape, axes) err = update_err(err, "cmaxf", _l2error(a.astype(np.complex64), b))
b=irfftn(rfftn(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,lastsize=lastsize,nthreads=nthreads) b = irfftn(rfftn(a.real.astype(np.float32), axes=axes, nthreads=nthreads),
err = _l2error(a.real,b) axes=axes, inorm=2, lastsize=lastsize, nthreads=nthreads)
if err > fmaxerr: err = update_err(err, "rmaxf", _l2error(a.real.astype(np.float32), b))
fmaxerr = err b = pypocketfft.separable_hartley(
print("fmaxerr:", fmaxerr, shape, axes) pypocketfft.separable_hartley(a.real, axes=axes, nthreads=nthreads),
b=ifftn(fftn(a.astype(np.complex64),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) axes=axes, inorm=2, nthreads=nthreads)
err = _l2error(a.astype(np.complex64),b) err = update_err(err, "hmax", _l2error(a.real, b))
if err > cmaxerrf: b = pypocketfft.genuine_hartley(
cmaxerrf = err pypocketfft.genuine_hartley(a.real, axes=axes, nthreads=nthreads),
print("cmaxerrf:", cmaxerrf, shape, axes) axes=axes, inorm=2, nthreads=nthreads)
b=irfftn(rfftn(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,lastsize=lastsize,nthreads=nthreads) err = update_err(err, "hmax", _l2error(a.real, b))
err = _l2error(a.real.astype(np.float32),b) b = pypocketfft.separable_hartley(
if err > fmaxerrf: pypocketfft.separable_hartley(
fmaxerrf = err a.real.astype(np.float32), axes=axes, nthreads=nthreads),
print("fmaxerrf:", fmaxerrf, shape, axes) axes=axes, inorm=2, nthreads=nthreads)
b=pypocketfft.separable_hartley(pypocketfft.separable_hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b))
err = _l2error(a.real,b) b = pypocketfft.genuine_hartley(
if err > hmaxerr: pypocketfft.genuine_hartley(a.real.astype(np.float32), axes=axes,
hmaxerr = err nthreads=nthreads),
print("hmaxerr:", hmaxerr, shape, axes) axes=axes, inorm=2, nthreads=nthreads)
b=pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = update_err(err, "hmaxf", _l2error(a.real.astype(np.float32), b))
err = _l2error(a.real,b)
if err > hmaxerr:
hmaxerr = err
print("hmaxerr:", hmaxerr, shape, axes)
b=pypocketfft.separable_hartley(pypocketfft.separable_hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real.astype(np.float32),b)
if err > hmaxerrf:
hmaxerrf = err
print("hmaxerrf:", hmaxerrf, shape, axes)
b=pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real.astype(np.float32),b)
if err > hmaxerrf:
hmaxerrf = err
print("hmaxerrf:", hmaxerrf, shape, axes)
err = dict()
while True: while True:
test() test(err)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment