Mailman 3 January 2011 - NumPy-Discussion

Help in speeding up accumulation in a matrix
by Nicolas SCHEFFER Jan. 30, 2011

Jan. 30, 2011

Hi all, First email to the list for me, I just want to say how grateful I am to have python+numpy+ipython etc... for my day to day needs. Great combination of software. Anyway, I've been having this bottleneck in one my algorithms that has been bugging me for quite a while. The objective is to speed this part up. I've been doing tons of optimization and parallel processing around that piece of code to get a decent run time. The problem is easy. You want to accumulate in a matrix, a weighted sum of other matrices. Let's call this function scale_and_add: def scale_and_add_re(R,w,Ms): (nb_add,mdim,j)=np.shape(Ms) for i in range(nb_add): R+=w[i]*Ms[i] return R This 'for' loop bugs me since I know this will slow things down. But the dimension of these things are so large that any attempt to vectorize this is slower and takes too much memory. I typically work with 1000 weights and matrices, matrices of dimension of several hundred. My current config is: In [2]: %timeit scale_and_add_re(R,w,Ms) 1 loops, best of 3: 392 ms per loop In [3]: w.shape Out[3]: (2000,) In [4]: Ms.shape Out[4]: (2000, 200, 200) I'd like to be able to double these dimensions. For instance I could use broadcasting by using a dot product %timeit dot(Ms.T,w) 1 loops, best of 3: 1.77 s per loop But this is i) slower ii) takes too much memory (btw, I'd really need an inplace dot-product in numpy to avoid the copy in memory, like the blas call in scipy.linalg. But that's for an other thread...) The matrices are squared and symmetric. I should be able to get something out of this, but I never found anything related to this in Numpy. I also tried a Cython reimplementation %timeit scale_and_add_reg(L1,w,Ms) 1 loops, best of 3: 393 ms per loop It brought nothing in speed. Here's the code @cython.boundscheck(False) def scale_and_add_reg(np.ndarray[float, ndim=2] R, np.ndarray[float, ndim=1] w, np.ndarray[float, ndim=3] Ms): return _scale_and_add_reg(R,w,Ms) @cython.boundscheck(False) cdef int _scale_and_add_reg(np.ndarray[float, ndim=2] R, np.ndarray[float, ndim=1] w, np.ndarray[float, ndim=3] Ms): cdef unsigned int mdim cdef int nb_add (nb_add,mdim,j)=np.shape(Ms) cdef unsigned int i for i from 0 <= i < nb_add: R+=w[i]*Ms[i] #for j in range(mdim): # for k in range(mdim): # R[j][k]+=w[i]*Ms[i][j][k] return 0 So here's my question if someone has time to answer it: Did I try anything possible? Should I move along and deal with this bottleneck? Or is there something else I didn't think of? Thanks for reading, keep up the great work! -n

6 7

[ANN] FEMTEC: Trac on open source scientific software
by Gael Varoquaux Jan. 30, 2011

Jan. 30, 2011

Hi list, This is just a note that an extra track at FEMTEC, a conference for computational methods in engineering and science, is open for open source scientific software. The organisers have a taste for Python, so if you want to submit a paper on numerical methods with Python, this is an excellent venue. Abstract submission is open till end of February. To submit you need to create an account and edit you profile. Gael ________________________________________________________________________________ The 3rd International Conference on Finite Element Methods in Engineering and Science (FEMTEC 2011, http://hpfem.org/events/femtec-2011/) will have a track on Open-source projects and Python in scientific computing. FEMTEC 2011 is co-organized by the University of Nevada (Reno), Desert Reseach Institute (Reno), Idaho National Laboratory (Idaho Falls, Idaho), and U.S. Army Engineer Research and Development Center (Vicksburg, Mississippi). The objective of the meeting is to strengthen the interaction between researchers who develop new computational methods, and scientists and engineers from various fields who employ numerical methods in their research. Specific focus areas of FEMTEC 2011 include, but are not limited to, the following: * Computational methods in hydrology, atmospheric modeling, and other earth sciences. * Computational methods in nuclear, mechanical, civil, electrical, and other engineering fields. * Mesh generation and scientific visualization. * Open-source projects and Python in scientific computing. Part of the conference will be a software afternoon featuring open source projects of participants. Proceedings Proceedings of FEMTEC 2011 will appear as a special issue of Journal of Computational and Applied Mathematics (2008 SCI impact factor 1.292), and additional high-impact international journals as needed.

1 0

Re: [Numpy-discussion] Wiener filter / decorrelation
by Sturla Molden Jan. 29, 2011

Jan. 29, 2011

Didn't I recommend the C++ and Fortran versions of 2nd edition? I particularly like the Fortran 90 edition as NumPy behaves like a vector machine as well. (The algorithms are explained in the Fortran 77 text, so they must be read together.) I'd like to warn against NR in C (all algorithms contain illegal C) and the third edition of NR (bloated and unorganised text, obfuscated C++), as well as known NR problems (e.g. bad PRNGs, unstable SVD code, slow FFTs, license). Another book I recommend is Golub and van Loan's text on matrix computations (1996). It explains the algorithms with array syntax (Matlab inspired), as well as citing LAPACK and BLAS routines relevant for the text. Sturla Den 28.01.2011 15:50, skrev denis: > Sturla, > what books do you like then, on any of NR's topics ? > Examples convince better than attacks. > Sure NR has its faults, but I like its style: de gustabus > cheers > -- denis > > On Jan 27, 11:47 pm, Sturla Molden<stu...(a)molden.no> wrote: > >> NR's third edition is utterly atrocious... >> extensive. It has many valuable details that should have been covered in >> previous versions, but it is presented in a way that makes me barf.

1 0

Should we make the master branch backward compatible.
by Charles R Harris Jan. 28, 2011

Jan. 28, 2011

Hi All, Mark Wiebe has proposed making the master branch backward compatible with 1.5. The argument for doing this is that 1) removing the new bits for new releases is a chore as the refactor schedule slips and 2) the new ABI isn't settled and keeping the current code in won't help with the merge. Mark thinks it is possible to keep the datetime types along with the new half types while restoring backward compatibility, and if so we won't lose anything by making the change. I'm in favor of this change, but I may have overlooked something. Thoughts? Chuck

5 6

SWIG examples from the Cookbook
by kant erburiad Jan. 28, 2011

Jan. 28, 2011

Hi, I'd like to ask for your help regarding the use of SWIG with numpy. ** problem description ** While I can compile successfully the examples provided in ./numpy/doc/swig I can't compile the first example provided in the Cookbook. http://www.scipy.org/Cookbook/SWIG_NumPy_examples "A simple ARGOUT_ARRAY1 example" I have this error: ezrange_wrap.c: In function ‘SWIG_AsVal_long’: ezrange_wrap.c:3227: error: initializer element is not constant ** hypothesis and ugly workaround ** Comparing the examples from ./numpy/doc/swig and from the Cookbook, I noticed the main difference lies in the fact the first are in C++ while the latter in plain C. I "converted" the example from the Coobook to look like a c++ project. Essentially I renamed the .c file with a .cxx extension and I modified setup.py accordingly. This time it compiles successfully and the module is usable in python. ** questions ** 1/ What should I do to compile the Cookbook example with a C compiler? 2/ If it appears swig is now only compatible with the C++ compiler, do you have a practical workaround to propose? (without renaming the C files with a .cxx extension) ? ** configuration ** debian testing i686 swig version 1.3.40 python-numpy version 1.5.1 (I also tried 1.4.1) python version 2.6 I hope I explained my problem clearly enough. Tell me if you need more details. Regards, KB

1 0

Re: [Numpy-discussion] Error in tanh for large complex argument
by Mark Bakker Jan. 28, 2011

Jan. 28, 2011

Good point, so we need a better solution that fixes all cases >> I'll file a ticket. >> >> Incidentally, if tanh(z) is simply programmed as >> >> (1.0 - exp(-2.0*z)) / (1.0 + exp(-2.0*z)) >This will overflow as z -> -\infty. The solution is probably to use a >different expression for Re(z) < 0, and to check how other libraries do >this in case the above still misses something. > > Pauli

3 3

incorrect behavior when complex number with zero imaginary part is multiplied by inf
by Mark Bakker Jan. 28, 2011

Jan. 28, 2011

When I multiply a complex number with inf, I get inf + inf*j: In [17]: inf * (1+1j) Out[17]: (inf+inf*j) Even when the imaginary part is really small: In [18]: inf * (1+1e-100j) Out[18]: (inf+inf*j) Yet when the imaginary part is zero (and it really is a real number), the imaginary part is nan: In [19]: inf * (1+0j) Out[19]: (inf+nan*j) That is not correct. It should really given (inf+0*j). (I know where it comes from, inf*0 is not defined, but in this case it is, as 1+0j is really a real number and inf is by definition real as well). If there is consensus I can file a ticket. Mark

2 1

Re: [Numpy-discussion] Developer NumPy list versus User NumPy list
by James A. Bednar Jan. 28, 2011

Jan. 28, 2011

| Date: Thu, 27 Jan 2011 16:32:22 -0700 | From: Charles R Harris <charlesr.harris(a)gmail.com> | | On Thu, Jan 27, 2011 at 4:23 PM, Robert Kern <robert.kern(a)gmail.com> wrote: | | > On Thu, Jan 27, 2011 at 17:17, Travis Oliphant <teoliphant(a)gmail.com> wrote: | > | > > Hey all, | > > | > > What is the thought about having two separate NumPy lists (one | > > for development discussions and one for user discussions)? | > | > We've resisted it for years. I don't think the split has done | > scipy much good. But that may just be my perspective because I'm | > subscribed to both and filter them both to the same folder. | | Me too. I don't think there is so much traffic that a distinction | needs to be made. I'm subscribed to the numpy digest, and I have 8 digest emails from yesterday (27 January), i.e. one single day, sitting in my inbox. These 8 digests represent who knows how many separate emails. If that is not heavy traffic, I really wouldn't know what is! As someone who uses numpy heavily (I manage a large numpy-based software project) but is not a numpy developer, I would very much appreciate having a separate user list. I can't bring myself to unsubscribe from the current list, for fear of not noticing some important new features, related packages, or serious issues, but sorting out those things from the rest of the posts does take significant work. None of my actual developers subscribe any more, as they found the volume of posts overwhelming, so I've sacrificed myself so that I can try to notice anything important and bring it to their attention. Anything that would help that would be greatly appreciated! Jim -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

2 1

Re: [Numpy-discussion] Error in tanh for large complex argument
by Mark Bakker Jan. 28, 2011

Jan. 28, 2011

I'll file a ticket. Incidentally, if tanh(z) is simply programmed as (1.0 - exp(-2.0*z)) / (1.0 + exp(-2.0*z)) the problem is fixed. Thanks, Mark [clip] > > Not for large complex values: > > > > In [85]: tanh(1000+0j) > > Out[85]: (nan+nan*j) > > Yep, it's a bug. Care to file a ticket? > > The implementation is just sinh/cosh, which overflows. > The fix is to provide an asymptotic expansion (sgn Re z), > although around the imaginary axis the switch is perhaps > somewhat messy to handle. > >

2 1

Developer NumPy list versus User NumPy list
by Travis Oliphant Jan. 27, 2011

Jan. 27, 2011

Hey all, What is the thought about having two separate NumPy lists (one for development discussions and one for user discussions)? -Travis

11 10