« back to all changes in this revision
Viewing changes to test/test_fftw3.py
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Package Import Robot
- Author(s): Jerome Kieffer
- Date: 2011-12-11 14:44:24 UTC
- Revision ID: package-import@ubuntu.com-20111211144424-n2fel2ww3dlajmuf
Tags: upstream-0.2.2
Import upstream version 0.2.2
added
removed
test/test_fftw3.py
1
import sys
2
sys.path.insert(0,'../build/lib')
3
sys.path.insert(0,'../build/lib.linux-x86_64-2.6')
4
import unittest
5
import time
6
import fftw3
7
import fftw3f
8
import fftw3l
9
from numpy.fft import fft,ifft,fftshift
10
import numpy as np
11
import fftw3.lib
12
import fftw3l.lib
13
import fftw3f.lib
14
import fftw3.planning
15
import fftw3l.planning
16
import fftw3f.planning
17
import os
18
19
h = 0.01
20
epsilon = 1e-1
21
beta = 1
22
N = 512
23
libs = [fftw3, fftw3f, fftw3l]
24
25
_complex = ['complex','singlecomplex','longcomplex']
26
_float = ['double','single','longdouble']
27
28
29
def fftw_propagation_aligned(N,repeats, lib, dtype):
30
t = np.linspace(-5,5,N)
31
dt = t[1]-t[0]
32
f = np.linspace(-1/dt/2.,1/dt/2.,N)
33
f = fftshift(f)
34
t = fftshift(t)
35
farray = lib.create_aligned_array(f.shape,dtype=dtype)
36
tarray = lib.create_aligned_array(t.shape,dtype=dtype)
37
fftplan = lib.Plan(tarray,farray,'forward')
38
ifftplan = lib.Plan(farray,tarray,'backward')
39
farray[:] = 0
40
tarray[:] = 0
41
tarray += np.exp(-t**2/0.5)
42
dispersion = np.exp(-1.j*h*beta*f)
43
ti = time.time()
44
for i in xrange(repeats):
45
fftplan()
46
farray *= dispersion/N
47
ifftplan()
48
to = time.time()-ti
49
return fftshift(t),fftshift(tarray),fftshift(f),fftshift(farray), to
50
51
def fftw_propagation(N,repeats,lib, dtype):
52
t = np.linspace(-5,5,N)
53
dt = t[1]-t[0]
54
f = np.linspace(-1/dt/2.,1/dt/2.,N)
55
f = fftshift(f)
56
t = fftshift(t)
57
farray = np.zeros(f.shape,dtype=dtype)
58
tarray = np.zeros(t.shape,dtype=dtype)
59
fftplan = lib.Plan(tarray,farray,'forward')
60
ifftplan = lib.Plan(farray,tarray,'backward')
61
farray[:] = 0
62
tarray[:] = 0
63
tarray += np.exp(-t**2/0.5)
64
dispersion = np.exp(-1.j*h*beta*f)
65
ti = time.time()
66
for i in xrange(repeats):
67
fftplan()
68
farray *= dispersion/N
69
ifftplan()
70
to = time.time()-ti
71
return fftshift(t),fftshift(tarray),fftshift(f),fftshift(farray), to
72
73
74
def np_propagation(N,repeats,dtype):
75
t = np.linspace(-5,5,N)
76
dt = t[1]-t[0]
77
f = np.linspace(-1/dt/2.,1/dt/2.,N)
78
f = fftshift(f)
79
t = fftshift(t)
80
tarray = np.zeros(t.shape,dtype)
81
tarray += np.exp(-t**2/0.5)
82
farray = np.zeros(tarray.shape,dtype)
83
dispersion = np.zeros(tarray.shape,dtype)
84
dispersion += np.exp(-1.j*h*beta*f)
85
ti = time.time()
86
for i in xrange(repeats):
87
farray = fft(tarray)
88
farray *= dispersion
89
tarray = ifft(farray)
90
to = time.time()-ti
91
return fftshift(t),fftshift(tarray),fftshift(f),fftshift(farray),to
92
93
94
def scipy_propagation(N,repeats, dtype):
95
try:
96
from scipy.fftpack import fft as sfft, ifft as sifft
97
except:
98
return None, None, None, None, np.NaN
99
t = np.linspace(-5,5,N)
100
dt = t[1]-t[0]
101
f = np.linspace(-1/dt/2.,1/dt/2.,N)
102
f = fftshift(f)
103
t = fftshift(t)
104
tarray = np.zeros(t.shape,dtype)
105
tarray += np.exp(-t**2/0.5)
106
farray = np.zeros(tarray.shape,dtype)
107
dispersion = np.zeros(tarray.shape,dtype)
108
dispersion += np.exp(-1.j*h*beta*f)
109
ti = time.time()
110
for i in xrange(repeats):
111
farray = sfft(tarray)
112
farray *= dispersion
113
tarray = sifft(tarray)
114
to = time.time()-ti
115
return fftshift(t),fftshift(tarray),fftshift(f),fftshift(farray),to
116
117
118
class ProductTestCase(unittest.TestCase):
119
120
def testSelect(self):
121
for lib in libs:
122
for plan in lib.lib._typelist:
123
if len(plan[1])>2:
124
plantype,(intype, outtype,length) = plan
125
shape = np.random.randint(2,5,length)
126
inputa = np.zeros(shape=shape, dtype=intype)
127
outputa = np.zeros(shape=shape, dtype=outtype)
128
else:
129
plantype, (intype,outtype) = plan
130
shape = np.random.randint(2,5,np.random.randint(4,8))
131
length = len(shape)
132
inputa = np.zeros(shape=shape, dtype=intype)
133
outputa = np.zeros(shape=shape, dtype=outtype)
134
func, name, types = lib.planning.select(inputa,outputa)
135
self.failUnless(name == plantype, "%s: select returned a "\
136
"wrong type for input array type=%s, output "\
137
"array type=%s, and dimension = %d" \
138
%(lib, inputa.dtype, outputa.dtype, length))
139
self.failUnless(func is getattr(lib.lib.lib, plantype), "%s: "\
140
"wrong library function for type %s"\
141
%(lib,plantype))
142
143
def testWisdom(self):
144
i=0
145
for lib in libs:
146
lib.forget_wisdom()
147
inputa = lib.create_aligned_array(1024,np.typeDict[_complex[i]])
148
outputa = lib.create_aligned_array(1024,np.typeDict[_complex[i]])
149
plan = lib.Plan(inputa,outputa,flags=['patient'])
150
soriginal = lib.export_wisdom_to_string()
151
lib.import_wisdom_from_string(soriginal)
152
lib.export_wisdom_to_file('test.wisdom')
153
lib.forget_wisdom()
154
lib.import_wisdom_from_file('test.wisdom')
155
os.remove('test.wisdom')
156
del inputa
157
del outputa
158
i+=1
159
160
def testPropagation(self):
161
#can only test for fftw3 because longdouble and single are not both implemented for scipy.fft and numpy.fft
162
Ns = [2**i for i in range(10,15)]
163
repeats = 2000
164
epsilon = 1e-3
165
times = []
166
for Nn in Ns:
167
t,A,f,B, ti = fftw_propagation(Nn,repeats, fftw3, _complex[0])
168
ft,fA,ff,fB, fti = fftw_propagation_aligned(Nn,repeats, fftw3, _complex[0])
169
nt,nA, nf, nB, nti = np_propagation(Nn,repeats, _complex[0])
170
st,sA, sf, sB, sti = scipy_propagation(Nn,repeats, _complex[0])
171
times.append((ti, fti,nti, sti))
172
self.failUnless(sum(abs(A)**2-abs(nA)**2)< epsilon, "Propagation "\
173
"of fftw3 and numpy gives "\
174
"different results")
175
self.failUnless(sum(abs(fA)**2-abs(nA)**2)< epsilon, "Propagation "\
176
"of aligned fftw3 and numpy gives "\
177
"different results")
178
print "Benchmark: %s" %fftw3
179
print "N fftw3 fftw3_aligned numpy scipy"
180
for i in range(len(Ns)):
181
print "%5d %5.2f %5.2f %5.2f %5.2f" %(Ns[i],\
182
times[i][0],\
183
times[i][1], \
184
times[i][2], \
185
times[i][3])
186
187
def test2D(self):
188
try:
189
from pylab import imread
190
except:
191
return
192
im = imread('Fourier2.png')
193
im = im[:,:,1]
194
a = np.zeros(im.shape, dtype=im.dtype)
195
a[:] = im[:]
196
b = np.zeros((im.shape[0],im.shape[1]/2+1),dtype=np.typeDict['singlecomplex'])
197
p = fftw3f.Plan(a,b,'forward')
198
ip = fftw3f.Plan(b,a,'backward')
199
p()
200
b/=np.prod(a.shape)
201
ip()
202
self.failUnless(a.sum()-im.sum() < epsilon, "2D fft and ifft did not "\
203
"reproduce the same image")
204
205
206
207
if __name__ == '__main__': unittest.main()
Loggerhead is a web-based interface for Breezy
Version: 2.0.1