Mailman 3 August 2007 - SciPy-Dev

ExtGen - a Python Extension module Generation tool - something very new and fresh
by Pearu Peterson Aug. 4, 2007

Aug. 4, 2007

Hi, I have started to work on an ExtGen package that is a high-level tool for constructing Python extension modules. I'll need it in the G3 F2PY project but ExtGen might have more general interest and usage, hence so early announcement. For more information and hello example, see http://www.scipy.org/ExtGen Enjoy, Pearu

2 2

Build problems
by Charles R Harris Aug. 4, 2007

Aug. 4, 2007

Hi All, Anyone see this before? atlas_blas_info: libraries lapack,blas,cblas,atlas not found in /usr/local/atlas/lib/ libraries lapack,blas,cblas,atlas not found in /usr/local/lib libraries lapack,blas,cblas,atlas not found in /usr/lib/sse2 libraries lapack,blas,cblas,atlas not found in /usr/lib NOT AVAILABLE $[charris@localhost scipy]$ ls /usr/local/atlas/lib libatlas.a libcblas.a libf77blas.a liblapack.a libptcblas.a libptf77blas.a libatlas.so libcblas.so libf77blas.so liblapack.so libptcblas.so libptf77blas.so How come the libraries aren't found? Numpy builds and runs fine. This is scipy revision 3220. Chuck

2 1

DocstringStandards webpage
by dmitrey Aug. 4, 2007

Aug. 4, 2007

hi all, there is a doctest format mentioned here http://projects.scipy.org/scipy/numpy/wiki/DocstringStandards can you attach an example in the web page? Also, it would be nice if any real func from scipy or numpy have been provided (as an example) in the web page. Regards, D.

2 1

max and argmax
by dmitrey Aug. 3, 2007

Aug. 3, 2007

hi all, in MATLAB [val ind] = max(arr) returns both max and index of the max element. In Python I use ind = argmax(arr) val = arr(ind) is it possible to avoid 2 lines, i.e. somehow to replace by single one? D.

2 2

Re: [SciPy-dev] FFTW performances in scipy and numpy
by David Cournapeau Aug. 3, 2007

Aug. 3, 2007

Steven G. Johnson wrote: > On Thu, 2 Aug 2007, David Cournapeau wrote: >> This is the way matlab works, right ? If I understand correctly, >> wisdoms are a way to compute plans "offline". So for example, if you >> compute plans with FFTW_MEASURE | FFTW_UNALIGNED, for inplace >> transforms and a set of sizes, you can record it in a wisdom, and >> reload it such as later calls with FFTW_MEASURE | FFTW_UNALIGNED will >> be fast ? > > Yes, although a wisdom file can contain as many saved plans as you want. > >> Anyway, all this sounds like it should be solved by adding a better >> infrastructure the current wrappers (ala matlab). > > I know a little about how the Matlab usage of FFTW works, and they are > definitely not getting the full performance you would get by calling > FFTW yourself from C etc. So they are not necessarily the gold standard. I certainly do not consider matlab as the gold standard. That's more a reason why the current situation to have worse performances than matlab for fft in scipy with fftw3 is not good. But I will work on a better cache mechanism with a user interface in python (the complex implementation of fft with fftw3 does not use copy anymore for a few days now, by the way). > If you malloc and then immediately free, most of the time the malloc > implementation is just going to re-use the same memory and so you will > get the same pointer over and over. So it's not a good test of malloc > alignment. Ah, this was stupid indeed. I should have checked the addresses returned by malloc. But then, playing a bit with the test program, I found that using size above ~ 17000 double starts to make the ratio of aligned data decreasing. This decreasing does not happen if I force malloc to use srbk and not mmap for big sizes: this is consistent with the fact the threshold for 32 bits for mmapping areas is 128 kb in gnu libc. Basically, areas which are allocated through mmap seem to be never 16 bytes aligned ! This is starting to go way beyond my knowledge... I thought mmap were page aligned, which means 16 bytes alignment. Maybe malloc does not return the pointer it got from mmap directly, but a shifted version for some reasons ? Maybe my test is flawed again ? I pasted it just below. For example, if you have N = 65384, the ratio is more about 10 % than 50 %; if you force not using mmap (M_MMAP_MAX = 0), then it goes back to ~50 %. regards, David --------------------------- #include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <time.h> #include <malloc.h> #define NARRAY 10000 #define N (1 << 16) int main(void) { void *a[NARRAY]; uintptr_t p; int i, nalign = 0; int st; /* default value, at least on my 32 bits ubuntu with glibc: 128 kb */ st = mallopt(M_MMAP_THRESHOLD, 128 * 1024); if (st == 0) { fprintf(stderr, "changing malloc option failed\n"); } st = mallopt(M_MMAP_MAX, NARRAY); if (st == 0) { fprintf(stderr, "changing malloc option failed\n"); } srand(time(NULL)); for (i = 0; i < NARRAY; ++i) { a[i] = malloc((rand() % N + 2) * sizeof(double)); p = (uintptr_t) a[i]; if (p % 16 == 0) ++nalign; } printf("%d/%d = %g%% are 16-byte aligned\n", nalign, NARRAY, nalign * 100.0/NARRAY); for (i = 0; i < NARRAY; ++i) free(a[i]); return 0; }

1 0

Re: [SciPy-dev] FFTW performances in scipy and numpy
by David Cournapeau Aug. 2, 2007

Aug. 2, 2007

Steven G. Johnson wrote: > On Wed, 1 Aug 2007, David Cournapeau wrote: >> - making numpy data 16 bytes aligned. This one is a bit tricky. I >> won't bother you with the details, but generally, numpy data may not >> be even "word aligned". Since some archs require some kind of >> alignement, there are some mechanisms to get aligned buffers from >> unaligned buffers in numpy API; I was thinking about an additional >> flag "SIMD alignement", since this could be quite useful for many >> optimized libraries using SIMD. But maybe this does not make sense, I >> have not yet thought enough about it to propose anything concrete to >> the main numpy developers. > > One possibility might be to allocate *all* arrays over a certain size > using posix_memalign, to get 16-byte alignnment, instead of malloc or > whatever you are doing now, since if the array is big enough the > overhead of alignment is negligible. It depends on the importance of > FFTs, of course, although there may well be other things where 16-byte > alignment will make SIMD optimization easier. I don't know anything > about Python/NumPy's memory management, however, so I don't know > whether this suggestion is feasible. from my understanding of numpy (quite limited, I cannot emphasize this enough), the data buffer can easily be made 16 bytes aligned, since they are always allocated using a macro which is calling malloc, but could replaced by posix_memalign. But then, there is another problem: many numpy arrays are just than one data buffer, are not always contiguous, etc... To interface with most C libraries, this means that a copy has to be made somewhat if the arrays are not contiguous or aligned; a whole infrastructure exists in numpy to do exactly this, but I do not know it really well yet. Whether replacing the current allocation with posix_memalign would enable the current infrastructure to always return 16 bytes aligned buffers is unclear to me, I should discuss this with people more knowledgeable than me. > > On x86, malloc is guaranteed to be 8-byte aligned, so barring > conspiracy the probability of 16-byte alignment should be 50%. On > x86_64, on the other hand, I've read that malloc is guaranteed to be > 16-byte aligned (for allocations > 16 bytes) by the ABI, so you should > have no worries in this case as long as Python uses malloc or similar. Se below for my findings about this problem. >> - I have tried FFTW_UNALIGNED + FFTW_ESTIMATE plans; unfortunately, >> I found that the performances were worse than using FFTW_MEASURE + >> copy (the copies are done into aligned buffers). I have since >> discover that this may be due to the totally broken architecture of >> my main workstation (a Pentium four): on my recent macbook (On linux, >> 32 bits, CoreDuo2), using no copy with FFTW_UNALIGNED is much better. > > If you use FFTW_UNALIGNED, then essentially FFTW cannot use SIMD (with > some exceptions). Whether a copy will be worth it for performance in > this case will depend upon the problem size. For very large problems, > I'm guessing it's probably not worth it, since SIMD has less of an > impact there and cache has a bigger impact. > >> - The above problem is fixable if we add a mechanisme to choose >> plans (ESTIMATE vs MEASURE vs ... I found that for 1d cases at least, >> ESTIMATE vs MEASURE is what really count performance wise). > > One thing that might help is to build in wisdom at least for a few > small power-of-two sizes, created at build time. e.g. wisdom for > sizes 128, 256, 512, 1024, 2048, 4096, and 8192 is < 2kB. > > (We also need to work on the estimator heuristics...this has proven to > be a hard problem in general, unfortunately.) This is the way matlab works, right ? If I understand correctly, wisdoms are a way to compute plans "offline". So for example, if you compute plans with FFTW_MEASURE | FFTW_UNALIGNED, for inplace transforms and a set of sizes, you can record it in a wisdom, and reload it such as later calls with FFTW_MEASURE | FFTW_UNALIGNED will be fast ? Anyway, all this sounds like it should be solved by adding a better infrastructure the current wrappers (ala matlab). > Your claim that numpy arrays are almost never 16-byte aligned strikes > me as odd; if true, it means that NumPy is doing something terribly > weird in its allocation. > > I just did a quick test, and on i386 with glibc (Debian GNU/Linux), > malloc'ing 10000 arrays of double variables of random size, almost > exactly 50% are 16-byte aligned as expected. And on x86_64, 100% are > 16-byte aligned (because of the abovementioned ABI requirement). Well, then I must be doing something wrong in my tests (I attached the test program I use to test pointers returned by malloc). This always return 1 for buffers allocated with fftw_malloc, but always 0 for buffers allocated with malloc on my machine, which is quite similar to yours (ubuntu, x86, glibc). Not sure what I am missing (surely something obvious). regards, David /* * test.c: allocate buffers of random size, and count the ones which are 16 bytes aligned */ #include <stdlib.h> #include <stdio.h> #include <string.h> #include <stddef.h> #include <math.h> #include <time.h> #ifdef FORCE_ALIGNED_MALLOC #include <fftw3.h> #define MALLOC(size) fftw_malloc((size)) #define FREE(ptr) fftw_free((ptr)) #else #define MALLOC(size) malloc((size)) #define FREE(ptr) free((ptr)) #endif /* Return 1 if the pointer x is 16 bytes aligned */ #define CHECK_SIMD_ALIGNMENT(x) \ ( ((ptrdiff_t)(x) & 0xf) == 0) const size_t MAXMEM = 1 << 26; int allocmem(size_t n); int main(void) { int i; int niter = 100000; size_t s; int acc = 0; double mean = 0, min = 1e100, max = 0; double r; /* not really random, but enough for here */ srand48(time(0)); for(i = 0; i < niter; ++i) { /* forcing the size to be in [1, MAXMEM] */ r = drand48() * MAXMEM; s = floor(r); if (s < 1) { s = 1; } else if (s > MAXMEM) { s = MAXMEM; } /* computing max and min allocated size for summary */ mean += r; if (max < s) { max = s; } if (min > s) { min = s; } /* allocating */ acc += allocmem(s); } fprintf(stdout, "ratio %d / %d; average size is %f (m %f, M %f)\n", acc, niter, mean / niter, min, max); return 0; } /* * Alloc n bytes, and return 0 if the memory allocated with this size was 16 * bytes aligned. */ int allocmem(size_t n) { char* tmp; int isal = 0; tmp = MALLOC(n); isal = CHECK_SIMD_ALIGNMENT(tmp); /* Not sure this is useful, but this should guarantee that tmp is really * allocated, and that the compiler is not getting too smart */ tmp[0] = 0; tmp[n-1] = 0; FREE(tmp); return isal; }

1 0

GSoC schedule question (letter for my mentors & Matthieu Brucher)
by dmitrey Aug. 1, 2007

Aug. 1, 2007

hi all, this is a letter primarily for my mentors & Matthieu Brucher according to the schedule proposed by my mentors I should work on chapter 2 2. Make existing openopt code usable by adding docstrings, unit tests, and documented sample scripts. Docstrings must conform to the standard. http://projects.scipy.org/scipy/numpy/wiki/CodingStyleGuidelines Tests must conform to the TestingGuidelines http://projects.scipy.org/scipy/scipy/wiki/TestingGuidelines Make sure ralg and lincher are fully functional, tested, and documented, completing as much as possible of your projected work on ralg (i.e., constraints to nonsmooth ralg solver (c<0, h=0, lb, ub, A, Aeq), as well as their derivatives). Documentation must be such as to allow other developers to understand and maintain the code, should you cease maintenance. Sample scripts should ensure that users have informative (and documented!) examples to work with. (7 days) So, one of the main problems for lincher (along with extern QP solver from cvxopt) is: a good line-search minimizer, that takes into account slope angle, is absent (we had discussed the problem some weeks ago in mailing lists). In scipy.optimize there is only one appropriate func: line_search. However, there are some problems with line_search docstring: line_search(f, myfprime, xk, pk, gfk, old_fval, old_old_fval, args=(), c1=0.0001, c2=0.90000000000000002, amax=50) Find alpha that satisfies strong Wolfe conditions. Uses the line search algorithm to enforce strong Wolfe conditions Wright and Nocedal, 'Numerical Optimization', 1999, pg. 59-60 For the zoom phase it uses an algorithm by Outputs: (alpha0, gc, fc) So, as you see, params are not described; especially I'm interested in old_fval and old_old_fval. So, there was a letter from NLPy developer about an alternative (link to an article was provided), also, Matthieu proposed using one of his own solvers. I think the importance of the auxiliary solver is very high and should be done in 1st order. My opinion about next thing to do is using an appropriate solver from Matthieu's package. However, Matthieu's syntax differs too much from openopt one. I think first of all there should be an openopt binding to Matthieu's package, that will allow for oo users to use same syntax: prob = NLP(...) r = prob.solve() I would implement the binding by myself, but I miss well-described API documentation of Matthieu's code. First of all I'm interested 1) what solvers does it contain? 2) how xtol, funtol, contol etc can be passed to the solvers? then, secondary (it can wait, maybe default parameters would be enough for now) 3) which params can I modify (like type of step (Armijo|Powell|etc), etc) BTW ralg is also missing a good line-search optimizer. It requires the one that finds solutions with slope angle > pi/2. But it can wait, it has one quite good and problem with lincher is more actual. So I think the next GSoC schedule step should be connection of 1-2 Matthieu's solvers (that take into account slope angle, like strong Wolfe conditions do) to native openopt syntax. Are you agree? if yes, my question to Matthieu: can you provide a description of some appropriate solvers from your package? and/or an example of usage? Regards, D.

2 3

Save complex value ndarray
by Xuedong Zhang Aug. 1, 2007

Aug. 1, 2007

Hi, Thanks for all people made up this wonderful software. I am just switching from matlab and still need to save some complex data into mat file format. The version I have (scipy.io.mio.savemat function) only save matV4 file and don't support complex data. I wonder if there is any recent work on SVN or I have to work out my own solution??? Thanks Xuedong

2 1