[Numpy-discussion] numpy ndarray questions

[Jochen]

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