Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Simon Perkins
ducc
Commits
3439205f
Commit
3439205f
authored
Jun 03, 2020
by
Martin Reinecke
Browse files
start renaming
parent
5e630a15
Changes
4
Hide whitespace changes
Inline
Side-by-side
pypocketfft/demos/bench.py
View file @
3439205f
import
numpy
as
np
import
ducc_0_1.
pypocketfft
as
pypocket
fft
import
ducc_0_1.
fft
as
ducc
fft
from
time
import
time
import
matplotlib.pyplot
as
plt
...
...
@@ -41,13 +41,12 @@ def measure_fftw_np_interface(a, nrepeat, nthr):
return
tmin
,
b
def
measure_pypocketfft
(
a
,
nrepeat
,
nthr
):
import
ducc_0_1.pypocketfft
as
ppf
def
measure_duccfft
(
a
,
nrepeat
,
nthr
):
tmin
=
1e38
b
=
a
.
copy
()
for
i
in
range
(
nrepeat
):
t0
=
time
()
b
=
ppf
.
c2c
(
a
,
out
=
b
,
forward
=
True
,
nthreads
=
nthr
)
b
=
duccfft
.
c2c
(
a
,
out
=
b
,
forward
=
True
,
nthreads
=
nthr
)
t1
=
time
()
tmin
=
min
(
tmin
,
t1
-
t0
)
return
tmin
,
b
...
...
@@ -110,7 +109,7 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="",
for
n
in
range
(
ntry
):
shp
=
np
.
random
.
randint
(
nmax
//
3
,
nmax
+
1
,
ndim
)
if
nice_sizes
:
shp
=
np
.
array
([
pypocket
fft
.
good_size
(
sz
)
for
sz
in
shp
])
shp
=
np
.
array
([
ducc
fft
.
good_size
(
sz
)
for
sz
in
shp
])
print
(
" {0:4d}/{1}: shape={2} ..."
.
format
(
n
,
ntry
,
shp
),
end
=
" "
,
flush
=
True
)
a
=
(
np
.
random
.
rand
(
*
shp
)
-
0.5
+
1j
*
(
np
.
random
.
rand
(
*
shp
)
-
0.5
)).
astype
(
tp
)
output
=
[]
...
...
@@ -130,8 +129,8 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="",
plt
.
show
()
funcs
=
(
measure_
pypocket
fft
,
measure_fftw
)
ttl
=
"
pypocket
fft/FFTW()"
funcs
=
(
measure_
ducc
fft
,
measure_fftw
)
ttl
=
"
ducc
fft/FFTW()"
ntry
=
100
nthr
=
1
nice_sizes
=
True
...
...
pypocketfft/demos/stress.py
View file @
3439205f
import
numpy
as
np
import
ducc_0_1.
pypocketfft
as
pypocket
fft
import
ducc_0_1.
fft
as
fft
def
_l2error
(
a
,
b
,
axes
):
...
...
@@ -7,23 +7,23 @@ def _l2error(a, b, axes):
def
fftn
(
a
,
axes
=
None
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
return
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
def
ifftn
(
a
,
axes
=
None
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
False
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
return
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
False
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
def
rfftn
(
a
,
axes
=
None
,
inorm
=
0
,
nthreads
=
1
):
return
pypocket
fft
.
r2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
nthreads
=
nthreads
)
return
fft
.
r2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
nthreads
=
nthreads
)
def
irfftn
(
a
,
axes
=
None
,
lastsize
=
0
,
inorm
=
0
,
nthreads
=
1
):
return
pypocket
fft
.
c2r
(
a
,
axes
=
axes
,
lastsize
=
lastsize
,
forward
=
False
,
inorm
=
inorm
,
nthreads
=
nthreads
)
return
fft
.
c2r
(
a
,
axes
=
axes
,
lastsize
=
lastsize
,
forward
=
False
,
inorm
=
inorm
,
nthreads
=
nthreads
)
nthreads
=
0
...
...
@@ -69,71 +69,71 @@ def test(err):
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
),
shape
)
b
=
pypocket
fft
.
separable_hartley
(
pypocket
fft
.
separable_hartley
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
),
b
=
fft
.
separable_hartley
(
fft
.
separable_hartley
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
),
axes
=
axes
,
inorm
=
2
,
nthreads
=
nthreads
)
err
=
update_err
(
err
,
"hmax"
,
_l2error
(
a
.
real
,
b
,
axes
),
shape
)
b
=
pypocket
fft
.
genuine_hartley
(
pypocket
fft
.
genuine_hartley
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
),
b
=
fft
.
genuine_hartley
(
fft
.
genuine_hartley
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
),
axes
=
axes
,
inorm
=
2
,
nthreads
=
nthreads
)
err
=
update_err
(
err
,
"hmax"
,
_l2error
(
a
.
real
,
b
,
axes
),
shape
)
b
=
pypocket
fft
.
separable_hartley
(
pypocket
fft
.
separable_hartley
(
b
=
fft
.
separable_hartley
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
genuine_hartley
(
pypocket
fft
.
genuine_hartley
(
a
.
real
.
astype
(
np
.
float32
),
axes
=
axes
,
b
=
fft
.
genuine_hartley
(
fft
.
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
),
shape
)
if
all
(
a
.
shape
[
i
]
>
1
for
i
in
axes
):
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
1
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
1
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
2
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
2
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
4
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dct
(
pypocket
fft
.
dct
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
4
),
b
=
fft
.
dct
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
1
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
1
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
2
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
2
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
4
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
b
=
pypocket
fft
.
dst
(
pypocket
fft
.
dst
(
a_32
.
real
,
axes
=
axes
,
nthreads
=
nthreads
,
type
=
4
),
b
=
fft
.
dst
(
fft
.
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
),
shape
)
...
...
pypocketfft/pypocketfft.cc
View file @
3439205f
...
...
@@ -652,7 +652,7 @@ out : int
void
add_pypocketfft
(
py
::
module
&
msup
)
{
using
namespace
pybind11
::
literals
;
auto
m
=
msup
.
def_submodule
(
"
pypocket
fft"
);
auto
m
=
msup
.
def_submodule
(
"fft"
);
m
.
doc
()
=
pypocketfft_DS
;
m
.
def
(
"c2c"
,
c2c
,
c2c_DS
,
"a"
_a
,
"axes"
_a
=
None
,
"forward"
_a
=
true
,
"inorm"
_a
=
0
,
"out"
_a
=
None
,
"nthreads"
_a
=
1
);
...
...
pypocketfft/test/test_pypocketfft.py
View file @
3439205f
import
ducc_0_1.
pypocketfft
as
pypocket
fft
import
ducc_0_1.
fft
as
fft
# import pyfftw
import
numpy
as
np
import
pytest
...
...
@@ -25,33 +25,33 @@ def _assert_close(a, b, epsilon):
def
fftn
(
a
,
axes
=
None
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
return
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
def
ifftn
(
a
,
axes
=
None
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
False
,
inorm
=
inorm
,
return
fft
.
c2c
(
a
,
axes
=
axes
,
forward
=
False
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
def
rfftn
(
a
,
axes
=
None
,
inorm
=
0
,
nthreads
=
1
):
return
pypocket
fft
.
r2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
return
fft
.
r2c
(
a
,
axes
=
axes
,
forward
=
True
,
inorm
=
inorm
,
nthreads
=
nthreads
)
def
irfftn
(
a
,
axes
=
None
,
lastsize
=
0
,
inorm
=
0
,
nthreads
=
1
):
return
pypocket
fft
.
c2r
(
a
,
axes
=
axes
,
lastsize
=
lastsize
,
forward
=
False
,
return
fft
.
c2r
(
a
,
axes
=
axes
,
lastsize
=
lastsize
,
forward
=
False
,
inorm
=
inorm
,
nthreads
=
nthreads
)
def
rfft_scipy
(
a
,
axis
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
r2r_fftpack
(
a
,
axes
=
(
axis
,),
real2hermitian
=
True
,
return
fft
.
r2r_fftpack
(
a
,
axes
=
(
axis
,),
real2hermitian
=
True
,
forward
=
True
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
def
irfft_scipy
(
a
,
axis
,
inorm
=
0
,
out
=
None
,
nthreads
=
1
):
return
pypocket
fft
.
r2r_fftpack
(
a
,
axes
=
(
axis
,),
real2hermitian
=
False
,
return
fft
.
r2r_fftpack
(
a
,
axes
=
(
axis
,),
real2hermitian
=
False
,
forward
=
False
,
inorm
=
inorm
,
out
=
out
,
nthreads
=
nthreads
)
...
...
@@ -194,7 +194,7 @@ def test_identity_r2(shp):
def
test_genuine_hartley
(
shp
):
rng
=
np
.
random
.
default_rng
(
np
.
random
.
SeedSequence
(
42
))
a
=
rng
.
random
(
shp
)
-
0.5
v1
=
pypocket
fft
.
genuine_hartley
(
a
)
v1
=
fft
.
genuine_hartley
(
a
)
v2
=
fftn
(
a
.
astype
(
np
.
complex128
))
v2
=
v2
.
real
+
v2
.
imag
assert_
(
_l2error
(
v1
,
v2
)
<
1e-15
)
...
...
@@ -204,7 +204,7 @@ def test_genuine_hartley(shp):
def
test_hartley_identity
(
shp
):
rng
=
np
.
random
.
default_rng
(
np
.
random
.
SeedSequence
(
42
))
a
=
rng
.
random
(
shp
)
-
0.5
v1
=
pypocket
fft
.
separable_hartley
(
pypocket
fft
.
separable_hartley
(
a
))
/
a
.
size
v1
=
fft
.
separable_hartley
(
fft
.
separable_hartley
(
a
))
/
a
.
size
assert_
(
_l2error
(
a
,
v1
)
<
1e-15
)
...
...
@@ -212,11 +212,11 @@ def test_hartley_identity(shp):
def
test_genuine_hartley_identity
(
shp
):
rng
=
np
.
random
.
default_rng
(
np
.
random
.
SeedSequence
(
42
))
a
=
rng
.
random
(
shp
)
-
0.5
v1
=
pypocket
fft
.
genuine_hartley
(
pypocket
fft
.
genuine_hartley
(
a
))
/
a
.
size
v1
=
fft
.
genuine_hartley
(
fft
.
genuine_hartley
(
a
))
/
a
.
size
assert_
(
_l2error
(
a
,
v1
)
<
1e-15
)
v1
=
a
.
copy
()
assert_
(
pypocket
fft
.
genuine_hartley
(
pypocket
fft
.
genuine_hartley
(
v1
,
out
=
v1
),
inorm
=
2
,
out
=
v1
)
is
v1
)
assert_
(
fft
.
genuine_hartley
(
fft
.
genuine_hartley
(
v1
,
out
=
v1
),
inorm
=
2
,
out
=
v1
)
is
v1
)
assert_
(
_l2error
(
a
,
v1
)
<
1e-15
)
...
...
@@ -225,7 +225,7 @@ def test_genuine_hartley_identity(shp):
def
test_genuine_hartley_2D
(
shp
,
axes
):
rng
=
np
.
random
.
default_rng
(
np
.
random
.
SeedSequence
(
42
))
a
=
rng
.
random
(
shp
)
-
0.5
assert_
(
_l2error
(
pypocket
fft
.
genuine_hartley
(
pypocket
fft
.
genuine_hartley
(
assert_
(
_l2error
(
fft
.
genuine_hartley
(
fft
.
genuine_hartley
(
a
,
axes
=
axes
),
axes
=
axes
,
inorm
=
2
),
a
)
<
1e-15
)
...
...
@@ -240,5 +240,5 @@ def testdcst1D(len, inorm, type, dtype):
itp
=
(
0
,
1
,
3
,
2
,
4
)
itype
=
itp
[
type
]
if
type
!=
1
or
len
>
1
:
# there are no length-1 type 1 DCTs
_assert_close
(
a
,
pypocketfft
.
dct
(
pypocket
fft
.
dct
(
a
,
inorm
=
inorm
,
type
=
type
),
inorm
=
2
-
inorm
,
type
=
itype
),
eps
)
_assert_close
(
a
,
pypocketfft
.
dst
(
pypocket
fft
.
dst
(
a
,
inorm
=
inorm
,
type
=
type
),
inorm
=
2
-
inorm
,
type
=
itype
),
eps
)
_assert_close
(
a
,
fft
.
dct
(
fft
.
dct
(
a
,
inorm
=
inorm
,
type
=
type
),
inorm
=
2
-
inorm
,
type
=
itype
),
eps
)
_assert_close
(
a
,
fft
.
dst
(
fft
.
dst
(
a
,
inorm
=
inorm
,
type
=
type
),
inorm
=
2
-
inorm
,
type
=
itype
),
eps
)
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment