Commit 804b4c6c authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'master' of gitlab.mpcdf.mpg.de:mtr/pypocketfft

parents bc355549 701fa516
import numpy as np
import pypocketfft
import pyfftw
from time import time
import matplotlib.pyplot as plt
import math
......
This diff is collapsed.
......@@ -2,8 +2,8 @@ import numpy as np
import pypocketfft
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
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):
......@@ -52,85 +52,89 @@ def test(err):
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
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))
err = update_err(err, "s4maxf", _l2error(a_32.real, b, axes))
err = dict()
......
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