Cleaning and fixing fft in scipy ?
Hi, When fixing the fft problem with prime number in numpy, I came again across the fft implementation in scipy, and related problems (suboptimal solution for fftw, see eg scipy ticket #1). I would like to improve the situation: - First, I think the module needs serious cleaning, as for now, it is a bunch of C files with many #define all accross the code, making it difficult to track things. I propose to split the sources by implementation (fft_fftw3.c, fft_fft2, fft_mkl, etc...); this can be done without any consequence on the implementation. - Then, improving fft where it is needed. Does this sound ok with scipy developers ? I should come up with a patch for zfft within today; even if I try to keep changes minimal, I would need people willing to test (I can easily test on linux with fftw2 and 3, and if necessary with mkl on the same platform. I would prefer avoiding testing on windows myself, as I have only a Japanese windows, making things extremely painful for me on this already painful platform). cheers, David
David Cournapeau wrote:
Hi,
When fixing the fft problem with prime number in numpy, I came again across the fft implementation in scipy, and related problems (suboptimal solution for fftw, see eg scipy ticket #1). I would like to improve the situation: - First, I think the module needs serious cleaning, as for now, it is a bunch of C files with many #define all accross the code, making it difficult to track things. I propose to split the sources by implementation (fft_fftw3.c, fft_fft2, fft_mkl, etc...); this can be done without any consequence on the implementation. - Then, improving fft where it is needed.
Does this sound ok with scipy developers ?
Go for it. Thank you! -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
David Cournapeau wrote:
Hi,
When fixing the fft problem with prime number in numpy, I came again across the fft implementation in scipy, and related problems (suboptimal solution for fftw, see eg scipy ticket #1). I would like to improve the situation: - First, I think the module needs serious cleaning, as for now, it is a bunch of C files with many #define all accross the code, making it difficult to track things.
Could you be more clear about where the problem in your eyes lies? There are multiple sources for the fft (original fftpack files + interface files to other fft libraries if the user has those installed). The ifdefs that I could find are just in the interface files to the "other" fft libraries that a user might have installed. I don't care if that is redone so that the setup.py file just uses different sources as opposed to defining pre-processor variables, but you might check with Pearu since he is the author of those interfaces. But, I also don't really see the problem with the way it is done, now.
I propose to split the sources by implementation (fft_fftw3.c, fft_fft2, fft_mkl, etc...); this can be done without any consequence on the implementation. - Then, improving fft where it is needed.
Please indicate what improvements you will be making. The fft's are performed by external libraries, I'm hesitant to start altering what those libraries are doing without clear justification. If you want to improve the interface to the optional fftw and mkl libraries that is one thing, but just saying you are going to "fix the fft implementation" is not re-assuring. -Travis
Travis Oliphant wrote:
Could you be more clear about where the problem in your eyes lies? There are multiple sources for the fft (original fftpack files + interface files to other fft libraries if the user has those installed).
The ifdefs that I could find are just in the interface files to the "other" fft libraries that a user might have installed. I don't care if that is redone so that the setup.py file just uses different sources as opposed to defining pre-processor variables, but you might check with Pearu since he is the author of those interfaces.
But, I also don't really see the problem with the way it is done, now.
I don't mean to discourage effort here. I'm just hesitant to start changing things without some kind of feedback from the original author (Pearu in this case). I just want to make sure he is aware of and agrees with changes that are being made. I know Pearu spent a fair bit of time creating the current fft interfaces. They may not be perfect but they are quite flexible and allow people to use fft's from multiple libraries. Most of the fft behavior comes from those multiple libraries. So, I just want to be clear about which fft implementation and interface is being modified and exactly what the problems are. -Travis
Travis Oliphant wrote:
Could you be more clear about where the problem in your eyes lies? There are multiple sources for the fft (original fftpack files + interface files to other fft libraries if the user has those installed).
The ifdefs that I could find are just in the interface files to the "other" fft libraries that a user might have installed. I don't care if that is redone so that the setup.py file just uses different sources as opposed to defining pre-processor variables, but you might check with Pearu since he is the author of those interfaces.
But, I also don't really see the problem with the way it is done, now.
I don't mean to discourage effort here. I'm just hesitant to start changing things without some kind of feedback from the original author (Pearu in this case). I agree, and that's exactly why I send this email.
I just want to make sure he is aware of and agrees with changes that are being made. I know Pearu spent a fair bit of time creating the current fft interfaces. They may not be perfect but they are quite flexible and allow people to use fft's from multiple libraries. I don't intend to change anything to the API. The reasoning is that scipy implementation of fft is suboptimal when using fftw3: it was really slow for some time, and I applied a quick fix, which is neither efficient or elegant (see ticket1 of scipy trac: http://projects.scipy.org/scipy/scipy/ticket/1). This patch took me a good time to get it right because of the fft sources organization. I understand that style is a matter of preferences more than anything else, and I am certainly not a reference, but I really think that with
Travis Oliphant wrote: the kind of following code, it is really difficult to understand what's going on: int i; complex_double *ptr = inout; #ifndef WITH_MKL #if defined(WITH_FFTW) fftw_plan plan = NULL; #endif #endif #if defined WITH_MKL DFTI_DESCRIPTOR_HANDLE desc_handle; #else double* wsave = NULL; #endif #ifdef WITH_FFTWORK coef_dbl* coef = NULL; #endif #ifndef WITH_MKL #ifdef WITH_DJBFFT int j; complex_double *ptrc = NULL; unsigned int *f = NULL; #endif #endif #ifdef WITH_FFTWORK if (ispow2le2e30(n)) { i = get_cache_id_zfftwork(n); coef = caches_zfftwork[i].coef; } else #endif #ifndef WITH_MKL #ifdef WITH_DJBFFT switch (n) { case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:; case 512:;case 1024:;case 2048:;case 4096:;case 8192: i = get_cache_id_zdjbfft(n); f = caches_zdjbfft[i].f; ptrc = (complex_double*)caches_zdjbfft[i].ptr; } if (f==0) #endif #endif #ifdef WITH_MKL desc_handle = caches_zmklfft[get_cache_id_zmklfft(n)].desc_handle; #elif defined WITH_FFTW plan = caches_zfftw[get_cache_id_zfftw(n,direction)].plan; #else wsave = caches_zfftpack[get_cache_id_zfftpack(n)].wsave; #endif I think it is much more readable to have one file with (fftw3) """ complex_double *ptr = inout; fftw_complex *ptrm = NULL; fftw_plan plan = NULL; int i; plan = caches_zfftw[get_cache_id_zfftw(n, dir)].plan; """ One with (fftw case) """ int i; complex_double *ptr = inout; fftw_plan plan = NULL; plan = caches_zfftw[get_cache_id_zfftw(n, direction)].plan; """ Etc... I guess the #ifdef solution was OK with one ot two implementation, but now, I fail to see how anyone can understand what the function does, which variables it is using, etc... To improve the fftw3 implementation, I think we must have a better caching system, and I don't see how I could do that with the way the source is organized right now.
Most of the fft behavior comes from those multiple libraries. So, I just want to be clear about which fft implementation and interface is being modified and exactly what the problems are.
No interface would be changed; actually, I wouldn't be surprised if the compiled object code would be exactly the same with my suggestion for cleaning up (To be sure that I don't introduce bugs, I actually compare each new file with the preprocessed original file with gcc -E). David
On Sat, April 28, 2007 10:16 am, David Cournapeau wrote:
Travis Oliphant wrote: <snip>
I don't mean to discourage effort here. I'm just hesitant to start changing things without some kind of feedback from the original author (Pearu in this case). I agree, and that's exactly why I send this email.
I certainly agree that scipy.fftpack implementation can be improved (as any piece of software). The current implementation was an attempt to support a large number of different fft implementations (they all have their own pros-and-cons) that had (sometimes very) different APIs. And hence the complexity of the current implementation. If I would start rewriting scipy.fftpack then I would probably choose a different approach to the current one too. However, I doubt that any other implementation with the same goal would be easy to understand or read to non-authors due to a variety of APIs of different fft implementations. In summary, I am not against to rewriting scipy.fftpack provided that 1) it is carried out in scipy.sandbox until it becomes more-or-less equivalent to the current scipy.fftpack in terms of features, unit-testing, supported fft backends, and documentation. 2) there are performance tests demonstrating that the new implementation outperforms the old one at least with the following backends: Fortran fftpack (that is default for most users), fftw2, and fftw3. Finally, I would also like to note that the efficiency of scipy.fftpack cacheing should be tested against both short sequences (N=64,..,512) as well as for long sequences (N>1024) as the same scheme could be very dependent of the sequence size when regarding perfomance. Best regards, Pearu
On Sat, April 28, 2007 1:06 pm, Pearu Peterson wrote:
In summary, I am not against to rewriting scipy.fftpack provided that 1) it is carried out in scipy.sandbox until it becomes more-or-less equivalent to the current scipy.fftpack in terms of features, unit-testing, supported fft backends, and documentation.
I now reread your original proposal and I think that this condition can be relaxed provided that after applying any patches to scipy.fftpack all fftpack unittests pass ok. Pearu
On Sat, April 28, 2007 1:06 pm, Pearu Peterson wrote:
In summary, I am not against to rewriting scipy.fftpack provided that 1) it is carried out in scipy.sandbox until it becomes more-or-less equivalent to the current scipy.fftpack in terms of features, unit-testing, supported fft backends, and documentation.
I now reread your original proposal and I think that this condition can be relaxed provided that after applying any patches to scipy.fftpack all fftpack unittests pass ok. Thanks for your email, Pearu. I would certainly not submit a patch which do not pass the unittest; problem being I cannot test all implementations (I can test fftpack/fftw2/fftw3 on Linux easily; I cannot get djbfft to compile on my machine; if necessary, I can install
Pearu Peterson wrote: the MKL; I don't know whether OS and compiler matters a lot for fftpack, as I don't intend to touch the setup part). I think that I was not really clear in my former emails: I do not suggest a rewrite, but a 2 steps improvements: - one which does not change anything to the code except its organization (one file for one fft API instead of the #ifdef) - once the first step is done and considered acceptable by the scipy developers, I would like to improve fftw3, such as it is at least on par with fftw2: right now, for fftw3, the arrays are copied twice. I think fftw3 requires a more sophisticated caching scheme, because with fftw3 you cannot use the same plan when the input/output are changed, at least with the basic API. For example, Octave is using such a scheme. Other implementations would be untouched. Cheers, David
Pearu
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
On Sat, April 28, 2007 1:52 pm, David Cournapeau wrote:
I think that I was not really clear in my former emails: I do not suggest a rewrite, but a 2 steps improvements: - one which does not change anything to the code except its organization (one file for one fft API instead of the #ifdef) - once the first step is done and considered acceptable by the scipy developers, I would like to improve fftw3, such as it is at least on par with fftw2: right now, for fftw3, the arrays are copied twice. I think fftw3 requires a more sophisticated caching scheme, because with fftw3 you cannot use the same plan when the input/output are changed, at least with the basic API. For example, Octave is using such a scheme. Other implementations would be untouched.
I am ok with both steps. Let us know when the first step is complete so that we can test that the interface is working for different backends. Pearu
Pearu Peterson wrote:
I am ok with both steps. Let us know when the first step is complete so that we can test that the interface is working for different backends.
I have something almost ready for review (complex 1d only), but I have a couple of questions: - My understanding is that except for djbfft backend, all other backend implement fft of any size (eg different than 2^n); so for djbfft backend, for size different than 2^n, it uses a diffferent backend. For example, if FFTW (v2) is selected and djbfft detected, then the zfft function is effectively using 2 backend: djbfft and fftw2. Would it be a big problem if in djbfft case, the sizes != 2^n are always handled through a fixed backend (fftpack I guess) ? - there is one backend called fftwork. I would like to check wether my patch at least compile for all backend, but I don't know this backend, and didn't find anything on the internet ? cheers, David
On 02/05/07, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
- My understanding is that except for djbfft backend, all other backend implement fft of any size (eg different than 2^n); so for djbfft backend, for size different than 2^n, it uses a diffferent backend. For example, if FFTW (v2) is selected and djbfft detected, then the zfft function is effectively using 2 backend: djbfft and fftw2. Would it be a big problem if in djbfft case, the sizes != 2^n are always handled through a fixed backend (fftpack I guess) ?
Well, probably not, as people using djbfft probably know what they are doing. But It does seem a bit drastic to drop from the fastest backend to the slowest if say fftw or MKL is also installed. Especially as this would provide lower performance than the current implementation. But why is this necessary anyway? It is only decided at compile time right? Or are you going the whole way and compiling all available back ends and then letting people choose dynamically? That would be a nice feature. Cheers, John
John Travers wrote:
On 02/05/07, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
- My understanding is that except for djbfft backend, all other backend implement fft of any size (eg different than 2^n); so for djbfft backend, for size different than 2^n, it uses a diffferent backend. For example, if FFTW (v2) is selected and djbfft detected, then the zfft function is effectively using 2 backend: djbfft and fftw2. Would it be a big problem if in djbfft case, the sizes != 2^n are always handled through a fixed backend (fftpack I guess) ?
Well, probably not, as people using djbfft probably know what they are doing. But It does seem a bit drastic to drop from the fastest backend to the slowest if say fftw or MKL is also installed. Especially as this would provide lower performance than the current implementation. If fftw or mkl is installed, I am not sure djbfft would be used anyway. It may have been the fastest fft a few years ago, but I cannot repeat any of the results given on djbfft with recent CPU (and by recent, I mean a pentium 4, that is already a few years old) and compilers. More fundamentally, if you care about efficiency fpr fft, you use 2^N size.
But why is this necessary anyway? It is only decided at compile time right? The goal of the reorganization was: - to make easier to read existing code - to make it easier to add additional backends
It was a bit complicated to support point 2 with only C macro in a reliable way (because of djbfft exception, like the ones I assume you put for mkl/djbfft). But after having thought a bit, I took a totally different approach: I am now generating the code with a python script (see my other email), such as adding a backend just consists in writing your implementation in a C file, and regenerating the python script. David
On May 2, 2007, at 05:02 , John Travers wrote:
On 02/05/07, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
- My understanding is that except for djbfft backend, all other backend implement fft of any size (eg different than 2^n); so for djbfft backend, for size different than 2^n, it uses a diffferent backend. For example, if FFTW (v2) is selected and djbfft detected, then the zfft function is effectively using 2 backend: djbfft and fftw2. Would it be a big problem if in djbfft case, the sizes != 2^n are always handled through a fixed backend (fftpack I guess) ?
Well, probably not, as people using djbfft probably know what they are doing.
Not necessarily :-) It's listed as a possible dependency, so people (like me) just install it. -- |>|\/|< /------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
On 28/04/07, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Thanks for your email, Pearu. I would certainly not submit a patch which do not pass the unittest; problem being I cannot test all implementations (I can test fftpack/fftw2/fftw3 on Linux easily; I cannot get djbfft to compile on my machine; if necessary, I can install the MKL; I don't know whether OS and compiler matters a lot for fftpack, as I don't intend to touch the setup part).
I can easily test the MKL version. I wrote the original patch for the MKL support and agree some improvements to the source would be nice. In fact MKL is only currently supported for complex ffts so I might get round to extending that at some point after you have done your improvements. Cheers, John
participants (6)
-
David Cournapeau -
David M. Cooke -
John Travers -
Pearu Peterson -
Robert Kern -
Travis Oliphant