Mailman 3 August 2008 - SciPy-User

KDE 2D, problem on axis and pylab.scatter
by F. B. 28 Aug '08

28 Aug '08

Hello I have a problem with this code. ------------------------------------------------------------------- import numpy as np from numpy import array from numpy import log import scipy.stats as stats from matplotlib.pyplot import imshow import pylab as pl import copy f=open("datain.txt","r") xyval=f.readlines() xvect=[] yvect=[] for i in range(len(xyval)): if (i< (len(xyval)/2) ): yvect.append(float(xyval[i])) else: xvect.append( log( float(xyval[i]) ) ) rvs=np.r_[[xvect],[yvect]] #yvalues=array(yvect) #xvalues=array(xvect) #xmin=xvalues.min() #xmax=xvalues.max() #ymin=yvalues.min() #ymax=yvalues.max() #stepx=xmax/256 #stepy=ymax/256 print "xmin: "+str(xmin)+" xmax: "+str(xmax)+" stepx: "+str(stepx)+" ymin: "+str(ymin)+" ymax: "+str(ymax)+" stepy: "+str(stepy) #x_flat=np.arange(0,xmax,stepx) #y_flat=np.arange(0,ymax,stepy) x_flat=np.arange(-25,-15,0.1) y_flat=np.arange(0,900,1) kde = stats.kde.gaussian_kde(rvs) #rvs.T = [[....],[....]] x,y = np.meshgrid(x_flat,y_flat) grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1) z = kde(grid_coords.T) z = z.reshape(len(y_flat),len(x_flat)) pl.hold(True) contplot=pl.contourf(z,20) ax=copy.deepcopy(pl.axis()) pl.scatter(xvect,yvect,c='y') pl.axis(ax) print pl.axis() pl.show() ------------------------------------------------------------------- I need to create a kernel density distribution of my data but the result that I obtain is strange. From readlines I read the x values that are in a range from 10**-7 to 10**-12, but for them (as you can see) I use the logarithm and then they go from -25 to -15. My y values are in range 50 to 1000 more or less. There are some problems, first of all, the range of the values. Using x_flat=np.arange(-25,-15,0.1) I have a problem with the value of the X axis. In the graph that I obtain, the minimum value in x is zero, and the maximum one is one hundred, but it is supposed to be between -25 and -15. (you can see that the length of x_flat is 100). The other problem that I have is the scatter plot. I should plot my data over the KDE but I cannot do that and I don't know why. However the shape of the KDE is that I expect. If someone could gently help me I would really appreciate it. Thanks and best regards F. B. _________________________________________________________________ Windows Live Mail: il programma gratuito per gestire tutta la tua posta! http://get.live.com/wlmail/overview

1 0

dgemm overwrite_c bug?
by Dan Goodman 28 Aug '08

28 Aug '08

Hi all, I'm interesting in doing a matrix multiplication in-place, and because numpy's dot function doesn't do that I had to go looking. After much pain I finally found the linalg.blas etc. stuff in scipy. However, I think that the dgemm function is ignoring the overwrite_c parameter, e.g.: from numpy import array from scipy import linalg gemm, = linalg.blas.get_blas_funcs(['gemm']) x=array([[1.,2.],[3.,4.]]) y=array([[5.,6.],[7.,8.]]) c=array([[8.,9.],[10.,11.]]) print gemm(1.0,x,y,1.,c,0,0,overwrite_c=1) print print c gives the output [[ 27. 31.] [ 53. 61.]] [[ 8. 9.] [ 10. 11.]] but the two matrices ought to be the same... Or have I misunderstood something? There seems to be a lack of documentation for these functions, so that's quite possible. Is there any way of getting access to an in-place matrix multiplication like this? Many thanks, Dan Goodman

1 0

Revisiting Log Arrays
by T J 27 Aug '08

27 Aug '08

Hi, A while back there was a discussion on including support for working with log arrays: http://www.mail-archive.com/numpy-discussion@scipy.org/msg08840.html It looked like some major steps were taken to this goal: http://www.mail-archive.com/numpy-discussion@scipy.org/msg08982.html but unfortunately, I didn't follow up from there. 1) The end of the second discussion made it seem like scipy.special was a good place for this. Has this been committed? I'd love to just do a 'svn update'. 2) If it has not been committed yet, I would still like to try it. What would I need to do to get this working? In particular, it looks like the patch was for numpy rather than scipy.special. 3) I'm curious how this works. Suppose I have a log array and that I multiply it (elementwise) by a normal array. Typically, 0 log 0 is defined to be 0. Would this be handled properly? Or would I need to post-process the array (that would be slow...)? 4) Eventually, it would be really nice to have logdot and the other functions mentioned in the second discussion. I'll do whatever I can to help move this along, but getting this committed and out there seems like a good next step. Thanks for the hard work already, T J

1 0

IFFT when FT is known on log scale
by Anand Patil 27 Aug '08

27 Aug '08

Hi all, I've got a Fourier transform stored for frequency values that are evenly spaced on the log scale, and I need to inverse Fourier transform it. I'm wondering whether there's a simple trick that would let me use the inverse fast Fourier transform. I need the log scale because the tails of the FT are long but I need a detailed representation near the origin. The domain spans many decades, so I can't just linearize. Any tips? The package NFFT3, http://www-user.tu-chemnitz.de/~potts/nfft/looks great but it's a pretty heavy weapon for this little problem. Thanks, Anand

1 0

SciPy, MPI and OpenMP
by Lorenzo Isella 27 Aug '08

27 Aug '08

Dear All, I have recently attended a crash course on MPI and OpenMP. The examples always involved C or Fortran code. Now, I have a thought: if working on a single processor, I hardly need to use pure C or pure Fortran. I usually write a Fortran code for the bottlenecks and compile it with f2py to create a python module I then import. Hence two questions: (1) Can I do something similar with many processors? E.g.: write a Python code, embed some compiled Fortran code which is supposed to run on many processors, get the results and come back to Python. Python--->Fortran on many processors--->back to Python. (2)Is it also possible to directly parallelize a Python code? I heard about thread locking in Python. I did some online research, there seems to be a lot of projects trying to combine Python and MPI/OpenMP, but many look rather "experimental". In particular, of course, I would like to hear about SciPy and parallel computing. Many thanks Lorenzo

4 3

Conditionally adding items to a list
by Roger Herikstad 26 Aug '08

26 Aug '08

Hi all, I have a prolem that I was wondering if anyone could come up with a good solution for. I basically have to lists of number and I want to add elements from one list to the other as long as the difference between the added element and all elements already in that list exceeds a certain threshold. The code I came up with is map(times1.append,ifilter(lambda(x): numpy.abs(x-numpy.array(times1)).min()>1000, times2)) but it quickly slows down if times2 becomes sufficiently long. I need to be able to do this with lists of 100,000++ elements. Does anyone know of a quicker way of doing this? Thanks! ~ Roger

4 5

interpolation speed problem
by Tribulations Parallèles 25 Aug '08

25 Aug '08

Hi everybody, I try to perform quick interpolation with Python. So far I have compared two solutions. First, the scipy solution: ##### import numpy as N from scipy import interpolate from time import time number = 1000000 a=N.arange( 0, number, 0.1) b=N.arange( number, 2*number, 0.1) f = interpolate.interp1d( a, b ) t_initial = time() for i in range(0,50): print "foo=%.5f" % f( 49999.5 ), t_final = time() print "\nTotal time =", t_final-t_initial ##### Now, the pygsl version (the two lines "astype(numx.float_)" are very important for the speed, these lines have been given by Pierre from the pygsl mailing list): ##### import pygsl.interpolation from time import time import numpy numx = pygsl._numobj number = 1000000 a=numpy.arange( 0, number ) b=numpy.arange( number, 2*number ) a = a.astype(numx.float_) b = b.astype(numx.float_) c = pygsl.interpolation.linear( len(a) ) pygsl.interpolation.linear.init( c , a , b) t_initial = time() for i in range(0,50): print "foo=%.5f" % c.eval(49999.5), t_final = time() print "\nTotal time =", t_final-t_initial #### The pygsl version takes on my machine 0.0004 s The scipy version takes 0.005 s. Is there some means to improve the scipy version? Thanks Julien

1 0

[ Python(x,y) ] New release : 2.0.3
by Pierre Raybaut 23 Aug '08

23 Aug '08

Hi all, As you may already know, Python(x,y) is a free scientific-oriented Python Distribution based on Qt and Eclipse providing a self-consistent scientific development environment. Release 2.0.3 is now available on http://www.pythonxy.com. (Full Edition, Basic Edition, and update patch will be available on monday on http://code.google.com/p/pythonxy) Changes history 08 -23 -2008 - Version 2.0.3 : * Added: o GDCM 2.0.8 (thanks to Mathieu Malaterre) - Grassroots DiCoM is a C++ library for dealing with DICOM medical files o pyExcelerator 0.6.3 - Generating Excel 97+ files, importing Excel 95+ files, support for UNICODE in Excel files, using variety of formatting features and printing options, Excel files and OLE2 compound files dumper o EasyGUI 0.83 - EasyGUI is a tiny Python module for very simple, very easy GUI programming * Updated: o PyDev 1.3.20 o Console 2.0.140 o PyQt 4.4.3 o Qt Eclipse Integration 1.4.1.1 (Qt help update) o VTK 5.0.4 o SymPy 0.6.2 o Cython 0.9.8.1.1 * Corrected: o maplotlib 0.98.3: compatibility issue with PyQt 4.4.x o Console 2: "Open console here..." now opens a command window instead of IPython-sh (the latter is less often used from Windows explorer, and can be opened afterwards thanks to Console2 multiple tabs management) Regards, Pierre Raybaut

1 0

Scilab to Scipy
by peter websdell 23 Aug '08

23 Aug '08

Hello all, I'm converting from using sclab/matlab to python. I'm finding it great so far, but I've discovered a problem that I don't understand. The following code displays a contour plot of the natural modes of a plate in scilab: ################### Lx=1; Ly=1; n=2; m=2; f=100; w=2*%pi*f; t=1; A=2; Kx=n*%pi/Lx; Ky=m*%pi/Ly; x=linspace(0,100); y=linspace(0,100); z=zeros(100,100); for i = 1:100 for j = 1:100 z(i,j) = A * sin(Kx*x(i)) * sin(Ky*y(j)) * %e^(%i*w*t); end end contour(x,y,z,20) ################## Now here's how I've done it in python: ################## from pylab import * Lx=1 Ly=1 n=2 m=2 f=100 w=2*pi*f t=1 A=2 Kx=n*pi/Lx Ky=m*pi/Ly x,y =mgrid[0:100,0:100] z=empty((100,100)) z=A * sin(Kx*x) * sin(Ky*y) * e**(1j*w*t) contour(x,y,z) show() ################### The result does plot a contour, but it is garbage. I have tried replecating the silly for loop approach, and also using the real of abs value of the result, but it is still garbage. Can anyone offer some advice as to why the two scripts produce different results? Thanks a lot, Pete.

3 8

How to use where
by Alexander Borghgraef 22 Aug '08

22 Aug '08

Hi all, I'm trying to figure out here how 'where' works exactly. I'm working on a list of vectors which I represent as a 2D array, and I'd like to remove the vectors which are out of bounds. So after some experimenting I got to basically this: listofvectors = ... # shape is ( 100, 2 ) bound = array( [ xmax, ymax ] ) inside = where( all( listofvectors < bound, axis = 1 ) ) # inside is ( array[ 1, 2, 4, 10, ... ] ) listofvectors = listofvectors[ inside, : ] # shape is ( 1, 100, 2 ) Ok, so I'm almost there. The change of shape is annoying, since I'll be using listofvectors in a loop to sample from an image. Not that I can't work around this, but I'm looking for a cleaner way. Also, I don't really get why numpy is doing this, anyone care to explain? Thanks. -- Alex Borghgraef

5 6