Merge branch 'new_interface' into 'master'

New interface

See merge request !12

bench.py

+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)
+ttl = "pypocketfft/FFTW()"
+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")

bench_nd.py

@@ -5,7 +5,7 @@ import pypocketfft
 from time import time
 import matplotlib.pyplot as plt
 
-nthreads=0
+nthreads=1
 def _l2error(a,b):
     return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
 
@@ -23,10 +23,10 @@ def bench_nd_fftn(ndim, nmax, ntry, tp, nrepeat, filename=""):
         tmin_pp=1e38
         for i in range(nrepeat):
             t0=time()
-            b=pypocketfft.fftn(a,nthreads=nthreads)
+            b=pypocketfft.c2c(a,nthreads=nthreads, forward=True)
             t1=time()
             tmin_pp = min(tmin_pp,t1-t0)
-        a2=pypocketfft.ifftn(b,fct=1./a.size)
+        a2=pypocketfft.c2c(b,inorm=2, forward=False)
         assert(_l2error(a,a2)<(2.5e-15 if tp=='c16' else 6e-7))
         res.append(tmin_pp/tmin_np)
     plt.title("t(pypocketfft / numpy 1.17), {}D, {}, max_extent={}".format(ndim, str(tp), nmax))

setup.py

@@ -21,9 +21,9 @@ if sys.platform == 'darwin':
     extra_compile_args += ['--stdlib=libc++', '-mmacosx-version-min=10.9']
     vars = distutils.sysconfig.get_config_vars()
     vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '')
-    python_module_link_args+=['-bundle']
+    python_module_link_args += ['-bundle']
 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']
 
 # if you don't want debugging info, add "-s" to python_module_link_args

stress.py

 import numpy as np
 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))
 
-nthreads=0
-cmaxerr=0.
-fmaxerr=0.
-cmaxerrf=0.
-fmaxerrf=0.
-hmaxerr=0.
-hmaxerrf=0.
-def test():
-    global cmaxerr, fmaxerr, hmaxerr, cmaxerrf, fmaxerrf, hmaxerrf
-    ndim=np.random.randint(1,5)
-    axlen=int((2**20)**(1./ndim))
-    shape = np.random.randint(1,axlen,ndim)
+
+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 not in err:
+        err[name] = value
+        print("{}: {}".format(name, value))
+    else:
+        if value > err[name]:
+            err[name] = value
+            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)
     np.random.shuffle(axes)
-    nax = np.random.randint(1,ndim+1)
+    nax = np.random.randint(1, ndim+1)
     axes = axes[:nax]
     lastsize = shape[axes[-1]]
-    fct = 1./np.prod(np.take(shape, axes))
-    a=np.random.rand(*shape)-0.5 + 1j*np.random.rand(*shape)-0.5j
-    b=pypocketfft.ifftn(pypocketfft.fftn(a,axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads)
-    err = _l2error(a,b)
-    if err > cmaxerr:
-        cmaxerr = err
-        print("cmaxerr:", cmaxerr, shape, axes)
-    b=pypocketfft.irfftn(pypocketfft.rfftn(a.real,axes=axes,nthreads=nthreads),axes=axes,fct=fct,lastsize=lastsize,nthreads=nthreads)
-    err = _l2error(a.real,b)
-    if err > fmaxerr:
-        fmaxerr = err
-        print("fmaxerr:", fmaxerr, shape, axes)
-    b=pypocketfft.ifftn(pypocketfft.fftn(a.astype(np.complex64),axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads)
-    err = _l2error(a.astype(np.complex64),b)
-    if err > cmaxerrf:
-        cmaxerrf = err
-        print("cmaxerrf:", cmaxerrf, shape, axes)
-    b=pypocketfft.irfftn(pypocketfft.rfftn(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,fct=fct,lastsize=lastsize,nthreads=nthreads)
-    err = _l2error(a.real.astype(np.float32),b)
-    if err > fmaxerrf:
-        fmaxerrf = err
-        print("fmaxerrf:", fmaxerrf, shape, axes)
-    b=pypocketfft.hartley(pypocketfft.hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads)
-    err = _l2error(a.real,b)
-    if err > hmaxerr:
-        hmaxerr = err
-        print("hmaxerr:", hmaxerr, shape, axes)
-    b=pypocketfft.hartley(pypocketfft.hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads)
-    err = _l2error(a.real.astype(np.float32),b)
-    if err > hmaxerrf:
-        hmaxerrf = err
-        print("hmaxerrf:", hmaxerrf, shape, axes)
+    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)
+    err = update_err(err, "cmax", _l2error(a, b))
+    b = ifftn(fftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
+              nthreads=nthreads)
+    err = update_err(err, "cmax", _l2error(a.real, b))
+    b = fftn(ifftn(a.real, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
+             nthreads=nthreads)
+    err = update_err(err, "cmax", _l2error(a.real, b))
+    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))
+    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))
+    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))
+    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))
+    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))
+    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))
+    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))
+
 
+err = dict()
 while True:
-    test()
+    test(err)