fftw3f.planning

1 # This file is part of PyFFTW. 2 # 3 # Copyright (C) 2009 Jochen Schroeder 4 # Email: jschrod@berlios.de 5 # 6 # PyFFTW is free software: you can redistribute it and/or modify 7 # it under the terms of the GNU General Public License as published by 8 # the Free Software Foundation, either version 3 of the License, or 9 # (at your option) any later version. 10 # 11 # PyFFTW is distributed in the hope that it will be useful, 12 # but WITHOUT ANY WARRANTY; without even the implied warranty of 13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 # GNU General Public License for more details. 15 # 16 # You should have received a copy of the GNU General Public License 17 # along with PyFFTW. If not, see <http://www.gnu.org/licenses/>. 18 19 import numpy as np 20 from numpy import typeDict 21 from lib import lib, _typelist, PyFile_AsFile, PyBuffer_FromReadWriteMemory, lib_threads 22 from ctypes import byref, c_double 23 24 fftw_flags = {'measure':0, 25 'destroy input': 1, 26 'unaligned': 2, 27 'conserve memory':4, 28 'exhaustive':8, 29 'preserve input': 16, 30 'patient': 32, 31 'estimate': 64} 32 33 realfft_type = {'halfcomplex r2c':0, 34 'halfcomplex c2r':1, 35 'discrete hartley':2, 36 'realeven 00':3, 37 'realeven 01':4, 38 'realeven 10':5, 39 'realeven 11':6, 40 'realodd 00':7, 41 'realodd 01':8, 42 'realodd 10':9, 43 'realodd 11':10} 44 45 46 fft_direction = {'forward' : -1, 'backward': 1} 47

48 -def create_aligned_array(shape, dtype=typeDict['singlecomplex'], boundary=16):

49 """Create an array which is aligned in memory 50 51 Parameters: 52 shape -- the shape of the array 53 dtype -- the dtype of the array (default=typeDict['singlecomplex']) 54 boundary -- the byte boundary to align to (default=16) 55 """ 56 N = np.prod(shape)*np.array(1,dtype).nbytes 57 tmp = np.zeros(N+boundary, dtype=np.uint8) 58 address = tmp.__array_interface__['data'][0] 59 offset = (boundary - address % boundary) 60 return tmp [offset:offset + N].view(dtype=dtype).reshape(shape)

61 62

63 -def execute(plan):

64 """Execute fftw-plan, i.e. perform Fourier transform on the arrays given 65 when the plan was created""" 66 lib.fftwf_execute(plan)

67

68 -def guru_execute_dft(plan, inarray, outarray):

69 """Guru interface: perform Fourier transform on two arrays, 70 outarray=fft(inarray) using the given plan. Important: This function 71 does not perform any checks on the array shape and alignment for 72 performance reasons. It is therefore crucial to only provide arrays 73 with the same shape, dtype and alignment as the arrays used for planning, 74 failure to do so can lead to unexpected behaviour and even python 75 segfaulting. 76 """ 77 lib.fftwf_execute_dft(plan, inarray, outarray)

78

79 -def destroy_plan(plan):

80 """Delete the given plan""" 81 if isinstance(plan,Plan): 82 del plan 83 else: 84 lib.fftwf_destroy_plan(plan)

85

86 -def select(inarray,outarray):

87 """From a given input and output np array select the appropriate 88 fftw3 plan to create.""" 89 if inarray.shape != outarray.shape: 90 if inarray.dtype == outarray.dtype: 91 raise TypeError, "Input array and output array must have the same "\ 92 "shape if they have the same dtype" 93 elif inarray.dtype == typeDict['singlecomplex'] and outarray.dtype == typeDict['single']: 94 inshape = list(outarray.shape) 95 inshape[-1] = inshape[-1]/2 + 1 96 if inarray.shape != tuple(inshape): 97 raise TypeError, "For complex to real transforms the complex "\ 98 "array must be of shape (n1 x n2 x...x "\ 99 "(n-1)/2 +1" 100 elif inarray.dtype == typeDict['single'] and outarray.dtype == typeDict['singlecomplex']: 101 outshape = list(inarray.shape) 102 outshape[-1] = outshape[-1]/2 + 1 103 if outarray.shape != tuple(outshape): 104 raise TypeError, "For real to complex transforms the complex "\ 105 "array must be of shape (n1 x n2 x...x "\ 106 "(n-1)/2 +1" 107 if inarray.dtype != typeDict['single'] and inarray.dtype != typeDict['singlecomplex']: 108 raise TypeError, "Input array has to be either floating point or"\ 109 " complex" 110 elif outarray.dtype != typeDict['single'] and outarray.dtype != typeDict['singlecomplex']: 111 raise TypeError, "Output array has to be either floating point "\ 112 "or complex" 113 i = 0 114 while(i < len(_typelist)): 115 name, types = _typelist[i] 116 if inarray.dtype != types[0]: 117 i += 8 118 continue 119 elif outarray.dtype != types[1]: 120 i += 4 121 continue 122 elif i in [3,7,11,15]: 123 return getattr(lib, name), name, types 124 elif len(inarray.shape) != types[2]: 125 i += 1 126 continue 127 else: 128 return getattr(lib, name), name, types

129

130 -def _create_complex2real_plan(inarray, outarray, flags):

131 """Internal function to create complex fft plan given an input and output 132 np array and the direction and flags integers""" 133 func, name, types = select(inarray,outarray) 134 #this is necessary because the r2c and c2r transforms always use the 135 #shape of the larger array (the real one) 136 if np.prod(inarray.shape) < np.prod(outarray.shape): 137 shape = outarray.shape 138 else: 139 shape = inarray.shape 140 141 if len(types) < 3: 142 plan = func(len(shape), np.asarray(shape, dtype=int), 143 inarray, outarray, flags) 144 if plan is None: 145 raise Exception, "Error creating fftwf plan %s for the given "\ 146 "parameters" %name 147 else: 148 return plan, name 149 elif types[2] == 1: 150 plan = func(shape[0], inarray, outarray, flags) 151 if plan is None: 152 raise Exception, "Error creating fftwf plan %s for the given "\ 153 "parameters" %name 154 else: 155 return plan, name 156 elif types[2] == 2: 157 plan = func(shape[0], shape[1], inarray, outarray, flags) 158 if plan is None: 159 raise Exception, "Error creating fftwf plan %s for the given "\ 160 "parameters" %name 161 else: 162 return plan, name 163 elif types[2] == 3: 164 plan = func(shape[0], shape[1], shape[2],inarray, outarray, flags) 165 if plan is None: 166 raise Exception, "Error creating fftwf plan %s for the given "\ 167 "parameters" %name 168 else: 169 return plan, name 170 else: 171 raise ValueError, 'the dimensions are not correct'

172 173

174 -def _create_complex_plan(inarray, outarray, direction, flags):

175 """Internal function to create complex fft plan given an input and output 176 np array and the direction and flags integers""" 177 func, name, types = select(inarray,outarray) 178 #this is necessary because the r2c and c2r transforms always use the 179 #shape of the larger array (the real one) 180 if np.prod(inarray.shape) < np.prod(outarray.shape): 181 shape = outarray.shape 182 else: 183 shape = inarray.shape 184 185 if len(types) < 3: 186 plan = func(len(shape), 187 np.asarray(shape, dtype=int), 188 inarray, outarray, direction, flags) 189 if plan is None: 190 raise Exception, "Error creating fftwf plan %s for the given "\ 191 "parameters" %name 192 else: 193 return plan, name 194 elif types[2] == 1: 195 plan = func(shape[0], inarray, outarray, direction, flags) 196 if plan is None: 197 raise Exception, "Error creating fftwf plan %s for the given "\ 198 "parameters" %name 199 else: 200 return plan, name 201 elif types[2] == 2: 202 plan = func(shape[0], shape[1], inarray, outarray,\ 203 direction, flags) 204 if plan is None: 205 raise Exception, "Error creating fftwf plan %s for the given "\ 206 "parameters" %name 207 else: 208 return plan, name 209 elif types[2] == 3: 210 plan = func(shape[0], shape[1], shape[2],\ 211 inarray, outarray, direction, flags) 212 if plan is None: 213 raise Exception, "Error creating fftwf plan %s for the given "\ 214 "parameters" %name 215 else: 216 return plan, name 217 else: 218 raise ValueError, 'the dimensions are not correct'

219

220 -def _create_real_plan(inarray, outarray, realtype, flags):

221 """Internal function to create real fft plan given an input and output 222 np array and the realtype and flags integers""" 223 if realtype == None: 224 raise ValueError, "Two real input arrays but no realtype list given" 225 func, name, types = select(inarray,outarray) 226 if len(types) < 3: 227 plan = func(len(inarray.shape), np.asarray(inarray.shape,dtype=int),\ 228 inarray, outarray, np.asarray(realtype), flags) 229 if plan is None: 230 raise Exception, "Error creating fftwf plan %s for the given "\ 231 "parameters" %name 232 else: 233 return plan, name 234 elif types[2] == 1: 235 plan = func(inarray.shape[0], inarray, outarray, realtype[0], flags) 236 if plan is None: 237 raise Exception, "Error creating fftwf plan %s for the given "\ 238 "parameters" %name 239 else: 240 return plan, name 241 elif types[2] == 2: 242 plan = func(inarray.shape[0], inarray.shape[1], inarray, outarray,\ 243 realtype[0], realtype[1], flags) 244 if plan is None: 245 raise Exception, "Error creating fftwf plan %s for the given "\ 246 "parameters" %name 247 else: 248 return plan, name 249 elif types[2] == 3: 250 plan = func(inarray.shape[0], inarray.shape[1],inarray.shape[2], \ 251 inarray, outarray, realtype[0], realtype[1], \ 252 realtype[2], flags) 253 if plan is None: 254 raise Exception, "Error creating fftwf plan %s for the given "\ 255 "parameters" %name 256 else: 257 return plan, name 258 else: 259 raise ValueError, 'the dimensions are not correct'

260

261 -def _create_plan(inarray, outarray, direction='forward', flags=['estimate'], 262 realtypes=None, nthreads=1):

263 """Internal function to create a complex fft plan given an input and output 264 np array and the direction and flags integers""" 265 if lib_threads is not None: 266 lib_threads.fftwf_plan_with_nthreads(nthreads) 267 elif nthreads > 1: 268 raise ValueError, "Cannot use more than 1 thread for non-threaded fftwf: %s" % (nthreads) 269 if inarray.dtype == np.typeDict['singlecomplex'] and \ 270 outarray.dtype == np.typeDict['singlecomplex']: 271 return _create_complex_plan(inarray,outarray, fft_direction[direction], 272 _cal_flag_value(flags)) 273 elif inarray.dtype == np.typeDict['singlecomplex'] or \ 274 outarray.dtype == np.typeDict['singlecomplex']: 275 return _create_complex2real_plan(inarray,outarray, 276 _cal_flag_value(flags)) 277 elif inarray.dtype == np.typeDict['single'] and \ 278 outarray.dtype == np.typeDict['single']: 279 return _create_real_plan(inarray,outarray, \ 280 [realfft_type[r] for r in realtypes],\ 281 _cal_flag_value(flags)) 282 else: 283 raise TypeError, "The input or output array has a dtype which is not supported by fftwf: %r, %r"\ 284 % (inarray.dtype, outarray.dtype)

285

286 -def _cal_flag_value(flags):

287 """Calculate the integer flag value from a list of string flags""" 288 ret = 0 289 for f in flags: 290 ret += fftw_flags[f] 291 return ret

292

293 -def print_plan(plan):

294 """Print a nerd-readable version of the plan to stdout""" 295 lib.fftwf_print_plan(plan)

296

297 -def fprint_plan(plan, filename):

298 """Print a nerd-readable version of the plan to the given filename""" 299 fp = open(filename, 'w') 300 c_fp = PyFile_AsFile(fp) 301 lib.fftwf_fprint_plan(plan, c_fp) 302 fp.close()

303

304 -class Plan(object):

305 """Object representing a fftw plan used to execute Fourier transforms in 306 fftw 307 308 Attributes: 309 shape -- the shape of the input and output arrays, i.e. the FFT 310 flags -- a list of the fft flags used in the planning 311 direction -- the direction of the FFT 312 ndim -- the dimensionality of the FFT 313 inarray -- the input array 314 outarray -- the output array 315 """

316 - def __init__(self, inarray=None, outarray=None, direction='forward', 317 flags=['estimate'], realtypes=None, create_plan=True, 318 nthreads = 1):

319 """Initialize the fftw plan. 320 Parameters: 321 inarray -- array to be transformed (default=None) 322 outarray -- array to contain the Fourier transform 323 (default=None) 324 If one of the arrays is None, the fft is considered to be 325 an inplace transform. 326 327 direction -- direction of the Fourier transform, forward 328 or backward (default='forward') 329 flags -- list of fftw-flags to be used in planning 330 (default=['estimate']) 331 realtypes -- list of fft-types for real-to-real ffts, this 332 needs to be given if both input and output 333 arrays are real (default=None) 334 create_plan -- weather to actually create the plan (default=True) 335 nthreads -- number of threads to be used by the plan, 336 available only for threaded libraries (default=1) 337 """ 338 339 self.flags = flags 340 self.direction = direction 341 self.realtypes = realtypes 342 self.nthreads = nthreads 343 if create_plan: 344 if inarray is None and outarray is None: 345 raise 'Need at least one array to create the plan' 346 elif outarray is None: 347 self.create_plan(inarray,inarray) 348 elif inarray is None: 349 self.create_plan(outarray,outarray) 350 else: 351 self.create_plan(inarray,outarray)

352

353 - def __set_shape(self,shape):

354 if len(shape)==1: 355 self.ndim = 1 356 self.N = tuple(shape) 357 358 elif len(shape) > 1: 359 self.ndim = len(shape) 360 self.N = tuple(shape) 361 else: 362 raise ValueError, 'shape must be at least one dimensional'

363

364 - def __get_shape(self):

365 return self.N

366 shape = property(__get_shape, __set_shape) 367

368 - def create_plan(self, inarray, outarray):

369 """Create the actual fftw-plan from inarray and outarray""" 370 self.plan, self.type_plan = _create_plan(inarray,outarray, 371 direction=self.direction, 372 flags=self.flags, 373 realtypes=self.realtypes, 374 nthreads=self.nthreads) 375 self.shape = inarray.shape 376 self.inarray = inarray 377 self.outarray = outarray

378

379 - def _get_parameter(self):

380 return self.plan

381 _as_parameter_ = property(_get_parameter) 382

383 - def __call__(self):

384 """Perform the Fourier transform outarray = fft(inarray) for 385 the arrays given at plan creation""" 386 self.execute()

387

388 - def execute(self):

389 """Execute the fftw plan, i.e. perform the FFT outarray = fft(inarray) 390 for the arrays given at plan creation""" 391 execute(self)

392

393 - def __del__(self):

394 destroy_plan(self)

395

396 - def guru_execute_dft(self,inarray,outarray):

397 """Guru interface: perform Fourier transform on two given arrays, 398 outarray=fft(inarray). Important: This method does not perform any 399 checks on the array shape and alignment for performance reasons. It is 400 therefore crucial to only provide arrays with the same shape, dtype and 401 alignment as the arrays used for planning, failure to do so can lead to 402 unexpected behaviour and possibly python segfaulting 403 """ 404 guru_execute_dft(self,inarray,outarray)

405

406 - def get_flops(self):

407 """Return an exact count of the number of floating-point additions, 408 multiplications, and fused multiply-add operations involved in 409 the plan's execution. The total number of floating-point 410 operations (flops) is add + mul + 2*fma, or add + mul + fma if 411 the hardware supports fused multiply-add instructions 412 (although the number of FMA operations is only approximate 413 because of compiler voodoo). 414 """ 415 add = c_double(0) 416 mul = c_double(0) 417 fma = c_double(0) 418 lib.fftwf_flops(self, byref (add), byref (mul), byref (fma)) 419 return add.value, mul.value, fma.value

420

Source Code for Module fftw3f.planning