Re: [Numpy-discussion] numpy ndarray questions
On 1/27/2009 6:03 AM, Jochen wrote: BTW memmap arrays have the same problem if I create a memmap array and later do something like a=a+1 all later changes will not be written to the file. = is Python's rebinding operator. a = a + 1 rebinds a to a different object. As for ndarray's, I'd like to point out the difference between a[:] = a + 1 and a = a + 1 S.M. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On 1/27/2009 1:26 AM, Jochen wrote:
a = fftw3.AlignedArray(1024,complex)
a = a+1
= used this way is not assignment, it is name binding.
It is easy to use function's like fftw_malloc with NumPy:
import ctypes
import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc
fftw_malloc.argtypes = [ctypes.c_ulong,]
fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype):
d = dtype()
address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong
if (address = 0): raise MemoryError, 'fftw_malloc returned NULL'
class Dummy(object): pass
d = Dummy()
d.__array_interface__ = {
'data' = (address, False),
'typestr' : dtype.str,
'descr' : dtype.descr,
'shape' : shape,
'strides' : strides,
'version' : 3
}
return numpy.asarray(d)
If you have to check for a particular alignment before calling fftw,
that is trivial as well:
def is_aligned(array, boundary):
address = array.__array_interface__['data'][0]
return not(address % boundary)
there a way that I get a different object type? Or even better is there
a way to prevent operations like a=a+1 or make them automatically
in-place operations?
a = a + 1 # rebinds the name 'a' to another array
a[:] = a + 1 # fills a with the result of a + 1
This has to do with Python syntax, not NumPy per se. You cannot overload
the behaviour of Python's name binding operator (=).
Sturla Molden
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On 1/27/2009 12:37 PM, Sturla Molden wrote: address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong if (address = 0): if (address == ): raise MemoryError, 'fftw_malloc returned NULL' Sorry for the typo. S.M. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On 1/27/2009 12:37 PM, Sturla Molden wrote: It is easy to use function's like fftw_malloc with NumPy: Besides this, if I were to write a wrapper for FFTW in Python, I would consider wrapping FFTW's Fortran interface with f2py. It is probably safer, as well as faster, than using ctypes. It would also allow the FFTW library to be linked statically to the Python extension, avoiding DLL hell. S.M. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On 1/27/2009 12:37 PM, Sturla Molden wrote:
import ctypes
import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc
fftw_malloc.argtypes = [ctypes.c_ulong,]
fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype):
d = dtype()
address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong
if (address = 0): raise MemoryError, 'fftw_malloc returned NULL'
class Dummy(object): pass
d = Dummy()
d.__array_interface__ = {
'data' = (address, False),
'typestr' : dtype.str,
'descr' : dtype.descr,
'shape' : shape,
'strides' : strides,
'version' : 3
}
return numpy.asarray(d)
Or if you just want to use numpy, aligning to a 16 byte boundary
can be done like this:
def aligned_array(N, dtype):
d = dtype()
tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8)
address = tmp.__array_interface__['data'][0]
offset = (16 - address % 16) % 16
return tmp[offset:offset+N].view(dtype=dtype)
S.M.
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 12:37 +0100, Sturla Molden wrote:
On 1/27/2009 1:26 AM, Jochen wrote:
a = fftw3.AlignedArray(1024,complex)
a = a+1
= used this way is not assignment, it is name binding.
It is easy to use function's like fftw_malloc with NumPy:
import ctypes
import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc
fftw_malloc.argtypes = [ctypes.c_ulong,]
fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype):
d = dtype()
address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong
if (address = 0): raise MemoryError, 'fftw_malloc returned NULL'
class Dummy(object): pass
d = Dummy()
d.__array_interface__ = {
'data' = (address, False),
'typestr' : dtype.str,
'descr' : dtype.descr,
'shape' : shape,
'strides' : strides,
'version' : 3
}
return numpy.asarray(d)
I actually do it slightly different, because I want to free the memory
using fftw_free. So I'm subclassing ndarray and in the __new__ function
of the subclass I create a buffer object from fftw_malloc using
PyBuffer_FromReadWriteMemory and pass that to ndarray.__new__. In
__del__ I check if the .base is a buffer object and do an fftw_free.
If you have to check for a particular alignment before calling fftw,
that is trivial as well:
def is_aligned(array, boundary):
address = array.__array_interface__['data'][0]
return not(address % boundary)
Yes I knew that, btw do you know of any systems which need alignment to
boundaries other than 16byte for simd operations?
there a way that I get a different object type? Or even better is there
a way to prevent operations like a=a+1 or make them automatically
in-place operations?
a = a + 1 # rebinds the name 'a' to another array
a[:] = a + 1 # fills a with the result of a + 1
I knew about this one, this is what I have been doing.
This has to do with Python syntax, not NumPy per se. You cannot overload
the behaviour of Python's name binding operator (=).
Yes I actually thought about it some more yesterday when going home and
realised that this would really not be possible. I guess I just have to
document it clearly so that people don't run into it.
Cheers
Jochen
Sturla Molden
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 14:16 +0100, Sturla Molden wrote:
On 1/27/2009 12:37 PM, Sturla Molden wrote:
import ctypes
import numpy
fftw_malloc = ctypes.cdll.fftw.fftw_malloc
fftw_malloc.argtypes = [ctypes.c_ulong,]
fftw_malloc.restype = ctypes.c_ulong
def aligned_array(N, dtype):
d = dtype()
address = fftw_malloc(N * d.nbytes) # restype = ctypes.c_ulong
if (address = 0): raise MemoryError, 'fftw_malloc returned NULL'
class Dummy(object): pass
d = Dummy()
d.__array_interface__ = {
'data' = (address, False),
'typestr' : dtype.str,
'descr' : dtype.descr,
'shape' : shape,
'strides' : strides,
'version' : 3
}
return numpy.asarray(d)
Or if you just want to use numpy, aligning to a 16 byte boundary
can be done like this:
def aligned_array(N, dtype):
d = dtype()
tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8)
address = tmp.__array_interface__['data'][0]
offset = (16 - address % 16) % 16
return tmp[offset:offset+N].view(dtype=dtype)
S.M.
Ah, I didn't think about doing it in python, cool thanks.
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 14:16 +0100, Sturla Molden wrote: def aligned_array(N, dtype): d = dtype() tmp = numpy.array(N * d.nbytes + 16, dtype=numpy.uint8) address = tmp.__array_interface__['data'][0] offset = (16 - address % 16) % 16 return tmp[offset:offset+N].view(dtype=dtype) Ah, I didn't think about doing it in python, cool thanks. Doing it from Python means you don't have to worry about manually deallocating the array afterwards. It seems the issue of 16-byte alignment has to do with efficient data alignment for SIMD instructions (SSE, MMX, etc). So this is not just an FFTW issue. I would just put a check for 16-byte alignment in the wrapper, and raise an exception (e.g. MemoryError) if the array is not aligned properly. Raising an exception will inform the user of the problem. I would not attempt to make a local copy if the array is errorneously aligned. That is my personal preference. Sturla Molden ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. This might help some: http://www.scipy.org/Subclasses Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Mon, 2009-01-26 at 19:25 -0600, Ryan May wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. This might help some: http://www.scipy.org/Subclasses Ryan Thanks, I had read about __array_finalize__, but not about __array_wrap__. I'm not quite sure if I understand how I can use this to get around my problem, can you explain? Cheers Jochen ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David Sorry maybe I wasn't quite clear, what I mean by the fftw way is creating a plan and executing the plan, instead of doing x=fft(y). As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. BTW memmap arrays have the same problem if I create a memmap array and later do something like a=a+1 all later changes will not be written to the file. BTW I just answered to you in the other thread on scipy-users as well, I'll try to just keep the rest of my posts here so people don't have to look at two lists if they want to follow things. Cheers Jochen __ _ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
Jochen wrote: On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David Sorry maybe I wasn't quite clear, what I mean by the fftw way is creating a plan and executing the plan, instead of doing x=fft(y). I guess I was not very clear, because my suggestion has nothing to do with the above. As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO. Another solution would be to work on numpy itself, so that it use aligned buffers, but that's obviously much longer term - that's just one of those things on my TODO list that I never took the time to do, David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David Sorry maybe I wasn't quite clear, what I mean by the fftw way is creating a plan and executing the plan, instead of doing x=fft(y). I guess I was not very clear, because my suggestion has nothing to do with the above. As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO. Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft. (I believe that's what you did in the old fftw bindings in scipy?). So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs. However I've been thinking about ways to make the whole thing more intuitive. At least now I can compare how much overhead I'm adding if I start to do checks in order to make the handling easier :) Another solution would be to work on numpy itself, so that it use aligned buffers, but that's obviously much longer term - that's just one of those things on my TODO list that I never took the time to do, I agree that would be nice, I'd think a couple of operations would probably benefit from that. David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
Jochen wrote: On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David Sorry maybe I wasn't quite clear, what I mean by the fftw way is creating a plan and executing the plan, instead of doing x=fft(y). I guess I was not very clear, because my suggestion has nothing to do with the above. As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO. Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft. Yes, but I don't see the problem with keeping old plans/creating new ones. It was a pain in C++, because I had to handle many backends, not just fftw, but in your case, doing it in python is not very difficult I think. (I believe that's what you did in the old fftw bindings in scipy?). that's what was done in scipy.fftpack originally, but not by me. So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs. Caching plans will not affect performances; on the contrary, you will be able to use more aggressively optimized plans if you add the ability to save/load plans to the disk (using wisdow mechanism). Checking for aligned arrays is one line in python, and should not affect much performances, specially if you do fft on big arrays, David ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy ndarray questions
On Tue, 2009-01-27 at 14:46 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 13:54 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 13:28 +0900, David Cournapeau wrote: Jochen wrote: On Tue, 2009-01-27 at 12:49 +0900, David Cournapeau wrote: Jochen wrote: Hi all, I just wrote ctypes bindings to fftw3 (see http://projects.scipy.org/pipermail/scipy-user/2009-January/019557.html for the post to scipy). Now I have a couple of numpy related questions: In order to be able to use simd instructions I create an ndarray subclass, which uses fftw_malloc to allocate the memory and fftw_free to free the memory when the array is deleted. This works fine for inplace operations however if someone does something like this: a = fftw3.AlignedArray(1024,complex) a = a+1 a.ctypes.data points to a different memory location (this is actually an even bigger problem when executing fftw plans), however type(a) still gives me class 'fftw3.planning.AlignedArray'. I can't comment about subclassing ndarrays, but I can give you a hint about aligned allocator problem: you could maintain two list of cached plans, automatically detect whether your arrays are aligned or not, and use the appropriate list of plans; one list is for aligned arrays, one for unaligned. Before removing support for fftw, I played with some C++ code to do exactly that. You can tell fftw to create plans for unaligned arrays by using FFTW_UNALIGNED flag: http://projects.scipy.org/scipy/scipy/browser/branches/refactor_fft/scipy/fftpack/backends/fftw3/src/zfft.cxx cheers, David Hi David, I have actually kept more closely to the fftw way of doing things, i.e. I create a plan python object from two arrays, it also stores the two arrays to prevent someone to delete the original arrays and then executing the plan. I am not sure I follow you when you say the fftw way: you can avoid having to store the arrays altogether while still using fftw plans, there is nothing unfftw about using plans that way. I think trying to guarantee that your arrays data buffers won't change is more complicated. David Sorry maybe I wasn't quite clear, what I mean by the fftw way is creating a plan and executing the plan, instead of doing x=fft(y). I guess I was not very clear, because my suggestion has nothing to do with the above. As you say the problem comes when you execute a plan on arrays which don't exist anymore, which causes python to segfault (I'm talking about using fftw_execute() not fftw_execute_dft). So yes essentially my problem is trying to ensure that the buffer does not change. I agree that if you just use fftw_execute, you will have segfaults: I am just not sure that your problem is to ensure that the buffer does not change :) You can instead handle relatively easily the case where the buffers are not aligned, while still using fftw API. The fftw_execute is just not an appropriate API for python code, IMHO. Ah ok, I think I understand you now :). I agree that using fftw_execute is not the most pythonic way of doing things. However every other way I could think of, you would need to do a number of checks and possibly create new plans, as well as keeping old plans around when doing a fft. Yes, but I don't see the problem with keeping old plans/creating new ones. It was a pain in C++, because I had to handle many backends, not just fftw, but in your case, doing it in python is not very difficult I think. (I believe that's what you did in the old fftw bindings in scipy?). that's what was done in scipy.fftpack originally, but not by me. So I decided to create fftw bindings first for maximum performance (I know, just doing the whole calculation in c is probably more appropriate if you want maximum performance ;), it also fits my needs. Caching plans will not affect performances; on the contrary, you will be able to use more aggressively optimized plans if you add the ability to save/load plans to the disk (using wisdow mechanism). Checking for aligned arrays is one line in python, and should not affect much performances, specially if you do fft on big arrays, I agree checking for aligned arrays does not affect performance much, I guess the lookup mechanism for the plans is what would take up the most time IMO. I'll definitely look into this though :) Cheers Jochen ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
