Mailman 3 March 2010 - NumPy-Discussion

Py3k: making a py3k compat header available in installed numpy for scipy
by David Cournapeau 30 Mar '10

30 Mar '10

Hi, I have worked on porting scipy to py3k, and it is mostly working. One thing which would be useful is to install something similar to npy_3kcompat.h in numpy, so that every scipy extension could share the compat header. Is the current python 3 compatibility header usable "in the wild", or will it still significantly change (this requiring a different, more stable one) ? cheers, David

4 4

Dealing with roundoff error
by Mike Sarahan 29 Mar '10

29 Mar '10

Hi all, I have run into some roundoff problems trying to line up some experimental spectra. The x coordinates are given in intervals of 0.1 units. I read the data in from a text file using np.loadtxt(). I think Robert's post here explains why the problem exists: http://mail.scipy.org/pipermail/numpy-discussion/2007-June/028133.html However, even linspace shows roundoff error: a=np.linspace(0.0,10.0,endpoint=False) b=np.linspace(0.1,10.1,endpoint=False) np.sum(a[1:]==b[:-1]) # Gives me 72, no 100 What is the best way to deal with it? Multiply the intervals by 10, then convert them to ints? Thanks, Mike

4 4

numpy.trapz() doesn't respect subclass
by Ryan May 29 Mar '10

29 Mar '10

Hi, I found that trapz() doesn't work with subclasses: http://projects.scipy.org/numpy/ticket/1438 A simple patch (attached) to change asarray() to asanyarray() fixes the problem fine. Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma

4 13

Why this Difference in Importing NumPy 1.2 vs 1.4?
by Wayne Watson 29 Mar '10

29 Mar '10

2 4

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 42, Issue 85
by Kevin Dunn 28 Mar '10

28 Mar '10

> Date: Sun, 28 Mar 2010 00:24:01 +0000 > From: Andrea Gavana <andrea.gavana(a)gmail.com> > Subject: [Numpy-discussion] Interpolation question > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: > <d5ff27201003271724o6c82ec75v225d819c84140b46(a)mail.gmail.com> > Content-Type: text/plain; charset=ISO-8859-1 > > Hi All, > > I have an interpolation problem and I am having some difficulties > in tackling it. I hope I can explain myself clearly enough. > > Basically, I have a whole bunch of 3D fluid flow simulations (close to > 1000), and they are a result of different combinations of parameters. > I was planning to use the Radial Basis Functions in scipy, but for the > moment let's assume, to simplify things, that I am dealing only with > one parameter (x). In 1000 simulations, this parameter x has 1000 > values, obviously. The problem is, the outcome of every single > simulation is a vector of oil production over time (let's say 40 > values per simulation, one per year), and I would like to be able to > interpolate my x parameter (1000 values) against all the simulations > (1000x40) and get an approximating function that, given another x > parameter (of size 1x1) will give me back an interpolated production > profile (of size 1x40). Andrea, may I suggest a different approach to RBF's. Realize that your vector of 40 values for each row in y are not independent of each other (they will be correlated). First perform a principal component analysis on this 1000 x 40 matrix and reduce it down to a 1000 x A matrix, called your scores matrix, where A is the number of independent components. A is selected so that it adequately summarizes Y without over-fitting and you will find A << 40, maybe 2 or 3. There are tools, such as cross-validation, that do this well enough. Then you can relate your single column of X to these independent column in A using a tool such as least squares: one least squares model per column in the scores matrix. This works because each column in the score vector is independent (contains totally orthogonal information) to the others. But I would be surprised if this works well enough, unless A = 1. But it sounds like your don't just have a single column in you X-variables (you hinted that the single column was just for simplification). In that case, I would build a projection to latent structures model (PLS) model that builds a single latent-variable model that simultaneously models the X-matrix, the Y-matrix as well as providing the maximal covariance between these two matrices. > Something along these lines: > > import numpy as np > from scipy.interpolate import Rbf > > # x.shape = (1000, 1) > # y.shape = (1000, 40) > > rbf = Rbf(x, y) > > # New result with xi.shape = (1, 1) --> fi.shape = (1, 40) > fi = rbf(xi) > > > Does anyone have a suggestion on how I could implement this? Sorry if > it sounds confused... Please feel free to correct any wrong > assumptions I have made, or to propose other approaches if you think > RBFs are not suitable for this kind of problems. > > Thank you in advance for your suggestions. > > Andrea. > > "Imagination Is The Only Weapon In The War Against Reality." > http://xoomer.alice.it/infinity77/ > > ==> Never *EVER* use RemovalGroup for your house removal. You'll > regret it forever. > http://thedoomedcity.blogspot.com/2010/03/removal-group-nightmare.html <==

1 1

ufunc improvements [Was: Warnings in numpy.ma.test()]
by Darren Dale 28 Mar '10

28 Mar '10

On Wed, Mar 17, 2010 at 10:16 PM, Charles R Harris <charlesr.harris(a)gmail.com> wrote: > > > On Wed, Mar 17, 2010 at 7:39 PM, Darren Dale <dsdale24(a)gmail.com> wrote: >> >> On Wed, Mar 17, 2010 at 8:22 PM, Charles R Harris >> <charlesr.harris(a)gmail.com> wrote: >> > >> > >> > On Wed, Mar 17, 2010 at 5:26 PM, Darren Dale <dsdale24(a)gmail.com> wrote: >> >> >> >> On Wed, Mar 17, 2010 at 5:43 PM, Charles R Harris >> >> <charlesr.harris(a)gmail.com> wrote: >> >> > On Wed, Mar 17, 2010 at 3:13 PM, Darren Dale <dsdale24(a)gmail.com> >> >> > wrote: >> >> >> On Wed, Mar 17, 2010 at 4:48 PM, Pierre GM <pgmdevlist(a)gmail.com> >> >> >> wrote: >> >> >> > On Mar 17, 2010, at 8:19 AM, Darren Dale wrote: >> >> >> >> >> >> >> >> I started thinking about a third method called __input_prepare__ >> >> >> >> that >> >> >> >> would be called on the way into the ufunc, which would allow you >> >> >> >> to >> >> >> >> intercept the input and pass a somehow modified copy back to the >> >> >> >> ufunc. The total flow would be: >> >> >> >> >> >> >> >> 1) Call myufunc(x, y[, z]) >> >> >> >> 2) myufunc calls ?.__input_prepare__(myufunc, x, y), which >> >> >> >> returns >> >> >> >> x', >> >> >> >> y' (or simply passes through x,y by default) >> >> >> >> 3) myufunc creates the output array z (if not specified) and >> >> >> >> calls >> >> >> >> ?.__array_prepare__(z, (myufunc, x, y, ...)) >> >> >> >> 4) myufunc finally gets around to performing the calculation >> >> >> >> 5) myufunc calls ?.__array_wrap__(z, (myufunc, x, y, ...)) and >> >> >> >> returns >> >> >> >> the result to the caller >> >> >> >> >> >> >> >> Is this general enough for your use case? I haven't tried to >> >> >> >> think >> >> >> >> about how to change some global state at one point and change it >> >> >> >> back >> >> >> >> at another, that seems like a bad idea and difficult to support. >> >> >> > >> >> >> > >> >> >> > Sounds like a good plan. If we could find a way to merge the first >> >> >> > two >> >> >> > (__input_prepare__ and __array_prepare__), that'd be ideal. >> >> >> >> >> >> I think it is better to keep them separate, so we don't have one >> >> >> method that is trying to do too much. It would be easier to explain >> >> >> in >> >> >> the documentation. >> >> >> >> >> >> I may not have much time to look into this until after Monday. Is >> >> >> there a deadline we need to consider? >> >> >> >> >> > >> >> > I don't think this should go into 2.0, I think it needs more thought. >> >> >> >> Now that you mention it, I agree that it would be too rushed to try to >> >> get it in for 2.0. Concerning a later release, is there anything in >> >> particular that you think needs to be clarified or reconsidered? >> >> >> >> > And >> >> > 2.0 already has significant code churn. Is there any reason beyond a >> >> > big >> >> > hassle not to set/restore the error state around all the ufunc calls >> >> > in >> >> > ma? >> >> > Beyond that, the PEP that you pointed to looks interesting. Maybe >> >> > some >> >> > sort >> >> > of decorator around ufunc calls could also be made to work. >> >> >> >> I think the PEP is interesting, but it is languishing. There were some >> >> questions and criticisms on the mailing list that I do not think were >> >> satisfactorily addressed, and as far as I know the author of the PEP >> >> has not pursued the matter further. There was some interest on the >> >> python-dev mailing list in the numpy community's use case, but I think >> >> we need to consider what can be done now to meet the needs of ndarray >> >> subclasses. I don't see PEP 3124 happening in the near future. >> >> >> >> What I am proposing is a simple extension to our existing framework to >> >> let subclasses hook into ufuncs and customize their behavior based on >> >> the context of the operation (using the __array_priority__ of the >> >> inputs and/or outputs, and the identity of the ufunc). The steps I >> >> listed allow customization at the critical steps: prepare the input, >> >> prepare the output, populate the output (currently no proposal for >> >> customization here), and finalize the output. The only additional step >> >> proposed is to prepare the input. >> >> >> > >> > What bothers me here is the opposing desire to separate ufuncs from >> > their >> > ndarray dependency, having them operate on buffer objects instead. As I >> > see >> > it ufuncs would be split into layers, with a lower layer operating on >> > buffer >> > objects, and an upper layer tying them together with ndarrays where the >> > "business" logic -- kinds, casting, etc -- resides. It is in that upper >> > layer that what you are proposing would reside. Mind, I'm not sure that >> > having matrices and masked arrays subclassing ndarray was the way to go, >> > but >> > given that they do one possible solution is to dump the whole mess onto >> > the >> > subtype with the highest priority. That subtype would then be >> > responsible >> > for casts and all the other stuff needed for the call and wrapping the >> > result. There could be library routines to help with that. It seems to >> > me >> > that that would be the most general way to go. In that sense ndarrays >> > themselves would just be another subtype with especially low priority. >> >> I'm sorry, I didn't understand your point. What you described sounds >> identical to how things are currently done. What distinction are you >> making, aside from operating on the buffer object? How would adding a >> method to modify the input to a ufunc complicate the situation? >> > > Just *one* function to rule them all and on the subtype dump it. No > __array_wrap__, __input_prepare__, or __array_prepare__, just something like > __handle_ufunc__. So it is similar but perhaps more radical. I'm proposing > having the ufunc upper layer do nothing but decide which argument type will > do all the rest of the work, casting, calling the low level ufunc base, > providing buffers, wrapping, etc. Instead of pasting bits and pieces into > the existing framework I would like to lay out a line of attack that ends up > separating ufuncs into smaller pieces that provide low level routines that > work on strided memory while leaving policy implementation to the subtype. > There would need to be some default type (ndarray) when the functions are > called on nested lists and scalars and I'm not sure of the best way to > handle that. > > I'm just sort of thinking out loud, don't take it too seriously. Thanks for the clarification. I think I see how this could work: if ufuncs were callable instances of classes, __call__ would find the input with highest priority and pass itself and the input to that object's __handle_ufunc__. Now it is up to __handle_ufunc__ to determine whether and how to modify the input, call some method on the ufunc (like execute) to perform the buffer operation, then __handle_ufunc__ performs the cast, deals with metadata and returns the result. I skipped a step: initializing the output buffer. Would that be rolled into the ufunc execution, or should it be possible for __handle_ufunc__ to access the initialized buffer before execution occurs(__array_prepare__)? I think it is important to be able to perform the cast and calculate metadata before ufunc execution. If an error occurs, an exception can be raised before the ufunc operates on the arrays, which can modifies the data in place. Darren

1 1

f2py: "could not crack entity declaration"
by David Warde-Farley 27 Mar '10

27 Mar '10

I decided to give wrapping this code a try: http://morrislab.med.utoronto.ca/~dwf/GLMnet.f90 I'm afraid my Fortran skills are fairly limited, but I do know that gfortran compiles it fine. f2py run on this file produces lots of errors of the form, Reading fortran codes... Reading file 'GLMnet.f90' (format:fix) Line #263 in GLMnet.f90:" real x(no,ni),y(no),w(no),vp(ni),ca(nx,nlam) 353" updatevars: could not crack entity declaration "ca(nx,nlam)353". Ignoring. Line #264 in GLMnet.f90:" real ulam(nlam),a0(nlam),rsq(nlam),alm(nlam) 354" updatevars: could not crack entity declaration "alm(nlam)354". Ignoring. Line #265 in GLMnet.f90:" integer jd(*),ia(nx),nin(nlam) 355" updatevars: could not crack entity declaration "nin(nlam)355". Ignoring. Line #289 in GLMnet.f90:" real x(no,ni),y(no),w(no),vp(ni),ulam(nlam) 378" updatevars: could not crack entity declaration "ulam(nlam)378". Ignoring. Line #290 in GLMnet.f90:" real ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam) 379" updatevars: could not crack entity declaration "alm(nlam)379". Ignoring. Line #291 in GLMnet.f90:" integer jd(*),ia(nx),nin(nlam) 380" updatevars: could not crack entity declaration "nin(nlam)380". Ignoring. Line #306 in GLMnet.f90:" call chkvars(no,ni,x,ju) 392" analyzeline: No name/args pattern found for li Is it the numbers that it is objecting to (I'm assuming these are some sort of punchcard thing)? Do I need to modify the code in some way to make it f2py-friendly? Thanks, David

3 5

numpy.array(arr.flat) mutates arr if arr.flags.fortran: bug?
by Zachary Pincus 27 Mar '10

27 Mar '10

Hello, I assume it is a bug that calling numpy.array() on a flatiter of a fortran-strided array that owns its own data causes that array to be rearranged somehow? Not sure what happens with a fancier-strided array that also owns its own data (because I'm not sure how to create one of those in python). This is from the latest svn version (2.0.0.dev8302) but was also present in a previous version too. Zach In [9]: a = numpy.array([[1,2],[3,4]]).copy('F') In [10]: a Out[10]: array([[1, 2], [3, 4]]) In [11]: list(a.flat) Out[11]: [1, 2, 3, 4] In [12]: a # no problem Out[12]: array([[1, 2], [3, 4]]) In [13]: numpy.array(a.flat) Out[13]: array([1, 2, 3, 4]) In [14]: a # this ain't right! Out[14]: array([[1, 3], [2, 4]]) In [15]: a = numpy.array([[1,2],[3,4]]).copy('C') In [16]: numpy.array(a.flat) Out[16]: array([1, 2, 3, 4]) In [17]: a Out[17]: array([[1, 2], [3, 4]])

3 3

StringIO test failure with Python3.1.2
by Bruce Southey 27 Mar '10

27 Mar '10

Hi, Wow, this is really impressive! I installed the svn numpy version '2.0.0.dev8300' with the latest Python 3.1.2 and it works! All the tests pass except: test_utils.test_lookfor I am guessing that it is this line as the other io imports do not have the period. from .io import StringIO ====================================================================== ERROR: test_utils.test_lookfor ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python3.1/site-packages/nose/case.py", line 177, in runTest self.test(*self.arg) File "/usr/local/lib/python3.1/site-packages/numpy/lib/tests/test_utils.py", line 10, in test_lookfor import_modules=False) File "/usr/local/lib/python3.1/site-packages/numpy/lib/utils.py", line 751, in lookfor cache = _lookfor_generate_cache(module, import_modules, regenerate) File "/usr/local/lib/python3.1/site-packages/numpy/lib/utils.py", line 852, in _lookfor_generate_cache from .io import StringIO ImportError: cannot import name StringIO ---------------------------------------------------------------------- Ran 2898 tests in 24.646s FAILED (KNOWNFAIL=5, errors=1) <nose.result.TextTestResult run=2898 errors=1 failures=0> Bruce

7 16

SciPy 2010 Tutorials: brainstorming and call for proposals
by Brian Granger 27 Mar '10

27 Mar '10

Greetings everyone, This year, there will be two days of tutorials (June 28th and 29th) before the main SciPy 2010 conference. Each of the two tutorial tracks (intro, advanced) will have a 3-4 hour morning and afternoon session both days, for a total of 4 intro sessions and 4 advanced sessions. The main tutorial web page for SciPy 2010 is here: http://conference.scipy.org/scipy2010/tutorials.html We are currently in the process of planning the tutorial sessions. You can help us in two ways: Brainstorm/vote on potential tutorial topics ============================================ To help us plan the tutorials, we have setup a web site that allow everyone in the community to brainstorm and vote on tutorial ideas/topics. The website for brainstorming/voting is here: http://conference.scipy.org/scipy2010/tutorialsUV.html The tutorial committee will use this information to help select the tutorials. Please jump in and let us know what tutorial topics you would like to see. Tutorial proposal submissions ============================= We are now accepting tutorial proposals from individuals or teams that would like to present a tutorial. Tutorials should be focused on covering a well defined topic in a hands on manner. We want to see tutorial attendees coding! We are pleased to offer tutorial presenters stipends this year for the first time: * 1 Session: $1,000 (half day) * 2 Sessions: $1,500 (full day) Optionally, part of this stipend can be applied to the presenter's registration costs. To submit a tutorial proposal please submit the following materials to 2010tutorials(a)scipy.org by April 15: * A short bio of the presenter or team members. * Which track the tutorial would be in (intro or advanced). * A short description and/or outline of the tutorial content. * A list of Python packages that attendees will need to have installed to follow along.

1 0

NumPy-Discussion March 2010