Mailman 3 April 2006 - SciPy-Dev

Arbitrary n-dimensional rebin routine
by Angus McMorland 12 Apr '06

12 Apr '06

Hi all, I'm not sure of the protocol here, so I'll try the list first... Here's an n-dimensional rebinning routine I've written (based loosely on IDLs congrid function), and hope it might be of some use. Perhaps it can get bundled with Gary Ruben's rebin function, currently on the wiki (http://scipy.org/Cookbook/Rebinning). I've thrown a few simple tests at the routine so far, with good results, but it could probably do with more field testing. The routine is currently optimised to take as long and as much memory as possible, so if someone with more knowledge of these things than me is interested to refine it a bit, that'd be great. As a scipy beginner, I'm also really interested to get some feedback on coding best-practices that I should have used and haven't. Cheers, Angus. def congrid(a, nud, method='neighbour', centre=False, minusone=False): '''Arbitrary resampling of source array to new dimension sizes. Currently only supports maintaining the same number of dimensions. Also doesn''t work for 1-D arrays unless promoted to shape (x,1). Based loosely on IDL''s congrid routine, which apparently originally came from a VAX/VMS routine of the same name. I''m not completely sure of the validity of using parallel 1-D interpolations repeated along each axis in succession, but the results are visually compelling. method: neighbour - closest value from original data the ''kinds'' supported by scipy.interpolate.interp1d centre: True - interpolation points are at the centres of the bins False - points are at the front edge of the bin minusone: For example- inarray.shape = (i,j) & new dimensions = (x,y) False - inarray is resampled by factors of (i/x) * (j/y) True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1) This prevents extrapolation one element beyond bounds of input array. ''' if not a.dtype.type in [n.typeDict['Float32'], n.typeDict['Float64']]: print "Converting to float" a = a.astype('Float32') if minusone: m1 = 1. else: m1 = 0. if centre: ofs = 0.5 else: ofs = 0. old = n.asarray( a.shape ) ndims = len( old ) if len( nud ) != ndims: print "Congrid: dimensions error. This routine currently only support rebinning to the same number of dimensions." return None nudr = n.asarray( nud ).astype('Float32') dimlist = [] if method == 'neighbour': for i in range( ndims ): base = nvals(i, nudr) dimlist.append( (old[i] - m1) / (nudr[i] - m1) \ * (base + ofs) - ofs ) cd = n.array( dimlist ) cdr = cd.round().astype( 'UInt16' ) nua = a[list( cdr )] return nua elif method in ['nearest','linear','cubic','spline']: # calculate new dims for i in range( ndims ): base = n.arange( nudr[i] ) dimlist.append( (old[i] - m1) / (nudr[i] - m1) \ * (base + ofs) - ofs ) # specify old dims olddims = [n.arange(i).astype('Float32') for i in list( a.shape )] # first interpolation - for ndims = any mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method ) nua = mint( dimlist[-1] ) trorder = [ndims - 1] + range( ndims - 1 ) for i in range( ndims - 2, -1, -1 ): nua = nua.transpose( trorder ) mint = scipy.interpolate.interp1d( olddims[i], nua, kind=method ) nua = mint( dimlist[i] ) if ndims > 1: # need one more transpose to return to original dimensions nua = nua.transpose( trorder ) return nua else: print "Congrid error: Unrecognized interpolation type.\n", \ "This routine currently only supports \'nearest\',\'linear\',\'cubic\', and \'spline\'." return None -- Angus McMorland email a.mcmorland(a)auckland.ac.nz mobile +64-21-155-4906 PhD Student, Neurophysiology / Multiphoton & Confocal Imaging Physiology, University of Auckland phone +64-9-3737-599 x89707 Armourer, Auckland University Fencing Secretary, Fencing North Inc.

2 4

Statistical review month : weighted histogram and cumfreq
by David Huard 11 Apr '06

11 Apr '06

I recently had to compute a weighted cumulative frequency distribution so I modified the scipy.stats.histogram and scipy.stats.cumfreq fonctions. I added a key in both function call, namely weight=None, where the default is simply uniform weights. I wanted to ask if this change would be welcome before submitting the patch. My concern is that the change modifies the result returned by the function. Presently, the histogram and cumfreq functions return integers arrays, the number of items lying in certain intervals. When these items are weighted, an integer count doesn't make much sense, and I normalized the histogram and cumfreq results. In other words, the new histogram function returns a float array of the frequency, instead of a count. I feel that having a normalized output is more pratical, but it would ruin existing code. There is always the possibility of creating a whistogram and wcumfreq functions, but this is not a pretty solution. I'd like your feedback about that. Cheers, David

2 2

ANN: SciPy 2006 Conference
by Travis N. Vaught 11 Apr '06

11 Apr '06

Greetings, The *SciPy 2006 Conference* is scheduled for August 17-18, 2006 at CalTech. A tremendous amount of work has gone into SciPy and Numpy over the past few months, and the scientific python community around these and other tools has truly flourished[1]. The Scipy 2006 Conference is an excellent opportunity to exchange ideas, learn techniques, contribute code and affect the direction of scientific computing with Python. Conference details are at http://www.scipy.org/SciPy2006 Keynote ------- Python language author Guido van Rossum (!) has agreed to be the Keynote speaker at this year's Conference. http://www.python.org/~guido/ Registration: ------------- Registration is now open. You may register early online for $100.00 at http://www.enthought.com/scipy06. Registration includes breakfast and lunch Thursday & Friday and a very nice dinner Thursday night. After July 14, 2006, registration will cost $150.00. Call for Presenters ------------------- If you are interested in presenting at the conference, you may submit an abstract in Plain Text, PDF or MS Word formats to abstracts(a)scipy.org -- the deadline for abstract submission is July 7, 2006. Papers and/or presentation slides are acceptable and are due by August 4, 2006. Tutorial Sessions ----------------- Several people have expressed interest in attending a tutorial session. The Wednesday before the conference might be a good day for this. Please email the list if you have particular topics that you are interested in. Here's a preliminary list: - Migrating from Numeric or Numarray to Numpy - 2D Visualization with Python - 3D Visualization with Python - Introduction to Scientific Computing with Python - Building Scientific Simulation Applications - Traits/TraitsUI Please rate these and add others in a subsequent thread to the SciPy-user mailing list. Perhaps we can pick 4-6 top ideas and recruit speakers as demand dictates. The authoritative list will be tracked here: http://www.scipy.org/SciPy2006/TutorialSessions Coding Sprints -------------- If anyone would like to arrive earlier (Monday and Tuesday the 14th and 15th of August), we can borrow a room on the CalTech campus to sit and code against particular libraries or apps of interest. Please register your interest in these coding sprints on the SciPy-user mailing list as well. The authoritative list will be tracked here: http://www.scipy.org/SciPy2006/CodingSprints Mailing list address: scipy-user(a)scipy.org Mailing list archives: http://dir.gmane.org/gmane.comp.python.scientific.user Mailing list signup: http://www.scipy.net/mailman/listinfo/scipy-user [1] Some stats: NumPy has averaged over 16,000 downloads per month Sept. 05 to March 06. SciPy has averaged over 3,800 downloads per month in Feb. and March 06. (both scipy and numpy figures do not include the 2000 instances per month downloaded as part of the Python Enthought Edition Distribution for Windows.)

1 0

scipy.sparse + umfpack + system_info
by Robert Cimrman 10 Apr '06

10 Apr '06

Hi all, I have added a basic umfpack detection support to numpy.distutils.system_info() - see example site.cfg below (and replace <UMFPACKv4.4_root>...). The Umfpack wrapper still resides in the sandbox, but its setup.py now uses the detection stuff. Now I would like to move the wrapper somewhere under scipy.sparse and I hit a problem: umpfack.py depends on scipy.sparse, since it uses the CSR/CSC matrix formats, but also scipy.sparse uses umfpack.py for solving linear equations.if present. How to get out of this circular dependency problem? Do you have any suggestions on how to organize the files? r. -- site.cfg: [umfpack] library_dirs = <UMFPACKv4.4_root>/UMFPACK/Lib:<UMFPACKv4.4_root>/AMD/Lib include_dirs = <UMFPACKv4.4_root>/UMFPACK/Include umfpack_libs = umfpack, amd

4 32

interp1d, bounds_error & fill_value
by Jonathan Taylor 05 Apr '06

05 Apr '06

just a note to mention that it seems that the docstring on interp1d is not quite correct: it states that NaN is returned if bounds_error=0. this is not quite true, as you can set fill_value (not mentioned in docstring but evident in code). if x and y are both int arrays, it does not seem to return NaN... >>> import scipy.interpolate >>> import numpy as N >>> >>> x=N.arange(20) >>> f=scipy.interpolate.interp1d(x,x,bounds_error=0, fill_value=0.) >>> print f(30) [0] >>> >>> g=scipy.interpolate.interp1d(x,x,bounds_error=0) >>> print g(30), N.isnan(g(30)) [-2147483648] [False] >>> >>> X=N.arange(0,20.) >>> h=scipy.interpolate.interp1d(X,X,bounds_error=0) >>> print h(30) [ nan] >>> >>> X=N.arange(0,20.) >>> i=scipy.interpolate.interp1d(X,x,bounds_error=0) >>> print i(30) [ nan] >>> >>> X=N.arange(0,20.) >>> j=scipy.interpolate.interp1d(x,X,bounds_error=0) >>> print j(30) [ nan] >>> -- ------------------------------------------------------------------------ I'm part of the Team in Training: please support our efforts for the Leukemia and Lymphoma Society! http://www.active.com/donate/tntsvmb/tntsvmbJTaylor GO TEAM !!! ------------------------------------------------------------------------ Jonathan Taylor Tel: 650.723.9230 Dept. of Statistics Fax: 650.725.8977 Sequoia Hall, 137 www-stat.stanford.edu/~jtaylo 390 Serra Mall Stanford, CA 94305

1 0

Stineman interpolation
by John Hunter 04 Apr '06

04 Apr '06

A couple of matplotlib users, who know much more about these matters than I do, put together this interpolation routine which they feel is superior and better behaved than many out there. With minor tweaks this may be suitable for scipy.interpolate (it currently lives in matplotlib.mlab svn) def slopes(x,y): """ SLOPES calculate the slope y'(x) Given data vectors X and Y SLOPES calculates Y'(X), i.e the slope of a curve Y(X). The slope is estimated using the slope obtained from that of a parabola through any three consecutive points. This method should be superior to that described in the appendix of A CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russel W. Stineman (Creative Computing July 1980) in at least one aspect: Circles for interpolation demand a known aspect ratio between x- and y-values. For many functions, however, the abscissa are given in different dimensions, so an aspect ratio is completely arbitrary. The parabola method gives very similar results to the circle method for most regular cases but behaves much better in special cases Norbert Nemec, Institute of Theoretical Physics, University or Regensburg, April 2006 Norbert.Nemec at physik.uni-regensburg.de (inspired by a original implementation by Halldor Bjornsson, Icelandic Meteorological Office, March 2006 halldor at vedur.is) """ # Cast key variables as float. x=nx.asarray(x, nx.Float) y=nx.asarray(y, nx.Float) yp=nx.zeros(y.shape, nx.Float) dx=x[1:] - x[:-1] dy=y[1:] - y[:-1] dydx = dy/dx yp[1:-1] = (dydx[:-1] * dx[1:] + dydx[1:] * dx[:-1])/(dx[1:] + dx[:-1]) yp[0] = 2.0 * dy[0]/dx[0] - yp[1] yp[-1] = 2.0 * dy[-1]/dx[-1] - yp[-2] return yp def stineman_interp(xi,x,y,yp=None): """ STINEMAN_INTERP Well behaved data interpolation. Given data vectors X and Y, the slope vector YP and a new abscissa vector XI the function stineman_interp(xi,x,y,yp) uses Stineman interpolation to calculate a vector YI corresponding to XI. Here's an example that generates a coarse sine curve, then interpolates over a finer abscissa: x = linspace(0,2*pi,20); y = sin(x); yp = cos(x) xi = linspace(0,2*pi,40); yi = stineman_interp(xi,x,y,yp); plot(x,y,'o',xi,yi) The interpolation method is described in the article A CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell W. Stineman. The article appeared in the July 1980 issue of Creative computing with a note from the editor stating that while they were not an academic journal but once in a while something serious and original comes in adding that this was "apparently a real solution" to a well known problem. For yp=None, the routine automatically determines the slopes using the "slopes" routine. X is assumed to be sorted in increasing order For values xi[j] < x[0] or xi[j] > x[-1], the routine tries a extrapolation. The relevance of the data obtained from this, of course, questionable... original implementation by Halldor Bjornsson, Icelandic Meteorolocial Office, March 2006 halldor at vedur.is completely reworked and optimized for Python by Norbert Nemec, Institute of Theoretical Physics, University or Regensburg, April 2006 Norbert.Nemec at physik.uni-regensburg.de """ # Cast key variables as float. x=nx.asarray(x, nx.Float) y=nx.asarray(y, nx.Float) assert x.shape == y.shape N=len(y) if yp is None: yp = slopes(x,y) else: yp=nx.asarray(yp, nx.Float) xi=nx.asarray(xi, nx.Float) yi=nx.zeros(xi.shape, nx.Float) # calculate linear slopes dx = x[1:] - x[:-1] dy = y[1:] - y[:-1] s = dy/dx #note length of s is N-1 so last element is #N-2 # find the segment each xi is in # this line actually is the key to the efficiency of this implementation idx = nx.searchsorted(x[1:-1], xi) # now we have generally: x[idx[j]] <= xi[j] <= x[idx[j]+1] # except at the boundaries, where it may be that xi[j] < x[0] or xi[j] > x[-1] # the y-values that would come out from a linear interpolation: sidx = nx.take(s, idx) xidx = nx.take(x, idx) yidx = nx.take(y, idx) xidxp1 = nx.take(x, idx+1) yo = yidx + sidx * (xi - xidx) # the difference that comes when using the slopes given in yp dy1 = (nx.take(yp, idx)- sidx) * (xi - xidx) # using the yp slope of the left point dy2 = (nx.take(yp, idx+1)-sidx) * (xi - xidxp1) # using the yp slope of the right point dy1dy2 = dy1*dy2 # The following is optimized for Python. The solution actually # does more calculations than necessary but exploiting the power # of numpy, this is far more efficient than coding a loop by hand # in Python yi = yo + dy1dy2 * nx.choose(nx.array(nx.sign(dy1dy2), nx.Int32)+1, ((2*xi-xidx-xidxp1)/((dy1-dy2)*(xidxp1-xidx)), 0.0, 1/(dy1+dy2),)) return yi # Here is some example code from pylab import figure, show, nx, linspace, stineman_interp x = linspace(0,2*nx.pi,20); y = nx.sin(x); yp = None xi = linspace(x[0],x[-1],100); yi = stineman_interp(xi,x,y,yp); fig = figure() ax = fig.add_subplot(111) ax.plot(x,y,'ro',xi,yi,'-b.') show()

1 0

[PATCH] d1mach problem
by Neal Becker 04 Apr '06

04 Apr '06

I think these patches will fix the problem with d1mach miscompiling with gcc4.1:

2 2

setup/config problem with scipy svn
by John Hunter 03 Apr '06

03 Apr '06

I just tried updating and building scipy from svn (which I've built on this platform recently in the last week) and am now getting python/svn/scipy> uname -a Linux peds-pc311 2.6.10-5-386 #1 Thu Aug 18 22:23:56 UTC 2005 i686 GNU/Linux python/svn/scipy> svn up At revision 1808. python/svn/scipy> rm -rf build python/svn/scipy> python setup.py build /usr/lib/python2.4/distutils/dist.py:236: UserWarning: Unknown distribution option: 'configuration' warnings.warn(msg) running build running config_fc python/svn/scipy> and nothing happens -- no build, no compile, etc... JDH

3 4

Re: [SciPy-dev] Statistics review months
by Pierre GM 03 Apr '06

03 Apr '06

Robert, Excellent initative, thanks a lot ! Before getting too involved, I have a question: should the functions support MaskedArrays (when possible) ? I think about var/std (already available for MA in this patch, that should be checked before inclusion [http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/ma-200603…]), or median, in particular...

1 0

ImportError: ndimage does not yet support 64-bit systems
by Nils Wagner 03 Apr '06

03 Apr '06

>>> from scipy import * import linsolve.umfpack -> failed: No module named _umfpack Traceback (most recent call last): File "<stdin>", line 1, in ? File "/usr/lib64/python2.4/site-packages/scipy/ndimage/__init__.py", line 33, in ? raise ImportError, "ndimage does not yet support 64-bit systems" ImportError: ndimage does not yet support 64-bit systems I know that I can disable config.add_subpackage('ndimage') in ~/scipy/Lib/setup.py but is there another option to fix the problem on 64bit systems ? Nils

1 0