Mailman 3 October 2010 - NumPy-Discussion

The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
by Rick Muller 27 Oct '10

27 Oct '10

Help! I'm having a problem in searching through the *elements* if a 2d array. I have a loop over a numpy array: n,m = G.shape print n,m for i in xrange(n): for j in xrange(m): print type(G), type(G[i,j]), type(float(G[i,j])) g = float(abs(G[i,j])) if g < cut: print i,j,G[i,j] However, I'm getting the error message: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() despite the fact that I'm not computing a truth value of an array. I'd like to just compare to G[i,j] or abs(G[i,j]), which should be a *single* value, and *not* an array. I even call 'float' here to make sure that I cast this to a normal python float. But I still get this error message. At this point, I suspect that I'm doing something dumb, rather than discovering some unique bug in numpy. Can anyone help me out? -- Rick Muller rpmuller(a)gmail.com 505-750-7557

3 3

adding 'order=' keyword arg to ravel and others
by Mark Wiebe 27 Oct '10

27 Oct '10

Because writing "arr.ravel('F')" doesn't seem as descriptive as "arr.ravel(order='F')", I wrote this simple patch. I added mention of the order='A' parameter to a few places it is relevant as well. Here's the branch on github: http://github.com/m-paradox/numpy/compare/master...ravel_keyword_arg Could someone review it for me? It also fixes the following ticket: http://projects.scipy.org/numpy/ticket/1581 Thanks, Mark

2 2

short circuit != ?
by Neal Becker 27 Oct '10

27 Oct '10

Is there a way to get a short circuit != ? That is, compare 2 arrays but stop as soon as the first element comparison fails? I'm assuming that np.all (a != b) will _not_ do this, but will first compare all elements.

8 10

NumPy FAQ and Python 3 support
by Peter 27 Oct '10

27 Oct '10

Hi all, This FAQ (found via Google) is out of date regarding Python 3 support: http://new.scipy.org/faq.html#does-numpy-currently-work-with-python-3-x-wha… Peter

1 0

Median method for ndarray
by Gökhan Sever 27 Oct '10

27 Oct '10

Hello, Is there a particular reason why there isn't a median method for ndarray? If not so, could we get one? Probably will go somewhere in this file: http://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/methods… I find it more practical to do array.mean() rather than np.mean(array), likewise would be useful having a similar method for median calculation. Thanks. -- Gökhan

2 1

Re: [Numpy-discussion] problems with numdifftools
by Nicolai Heitz 27 Oct '10

27 Oct '10

Am 26.10.2010 12:38, schrieb josef.pktd(a)gmail.com: > On Tue, Oct 26, 2010 at 3:28 PM, Pauli Virtanen<pav(a)iki.fi> wrote: >> Tue, 26 Oct 2010 12:16:53 -0700, Nicolai Heitz wrote: >>> I am not sure if you are the right persons to contact but if not I would >>> appreciate a short notice and maybe an address where I can find help. I >>> already posted this this message in an other python mailing list and >>> they forwarded me to this list and told me that I might find help here. >> Some comments: >> >> 1) General advice: When doing numerics, it's generally a good idea to use >> units natural for the problem, and not SI ones. The numbers the computer >> sees should be of order 1. Not all numerical algorithms are scale >> invariant. > I also think it will be difficult to get good numbers by numerical > differenttiation with this scaling. > I knew before that the scaling might causes some problems but I didn't expect them to arise on that level. I thought numdifftools should be robust in those typical physical units/scales. I can transform it into natural units but this is not the best solution and I still hope for somebody to come up with a solution for SI units. >> 2) The scipy-user list might be even more appropriate: >> >> http://mail.scipy.org/mailman/listinfo/scipy-user I contacted them already but they didn't responded so far and I was forwarded to that list which was supposed to be more appropriated. >> 3) The numdifftools authors might be more knowledgeable about their >> software: >> >> http://code.google.com/p/numdifftools/ >> >> If you are sure it's a bug in numdifftools, click the issues tab and >> write a report. Be sure to include simple test case (like the one you >> attached). > Per might be reading this or scipy-user. > > My guess would be that it is not a bug but numerical precision > problems in a difficult use case. The question is however useful, > because I haven't seen much discussion yet about a robust use of > numdifftools. > > Josef > I am not sure if it is a bug either. I mean for most of the numbers I tested the code (not having a magnitude of e-28) like 2.5, 1, 1e-10 numdifftools works all fine. My question is more like 1) Can I make it run/fix it, so that it is also going to work for the SI scaling? 2) How can I be sure that increasing the number of ions or adding a somehow more complicated term to the potential energy is not causing the same problems even in natural units? 3) In which range is numdifftools working properly. Nicolai >> -- >> Pauli Virtanen >> >> _______________________________________________ >> NumPy-Discussion mailing list >> NumPy-Discussion(a)scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion(a)scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion

3 2

problems with numdifftools
by Nicolai Heitz 26 Oct '10

26 Oct '10

Hey, I am not sure if you are the right persons to contact but if not I would appreciate a short notice and maybe an address where I can find help. I already posted this this message in an other python mailing list and they forwarded me to this list and told me that I might find help here. I am using Python 2.6 and the following packages: 1) numpy 2) scipy 3) numdifftools I am a physicist and started to use Python 2-3 weeks ago. I want to use Python to calculate the eigenvalues of the Hamiltonian given in the code. Therefore I have to do the following steps: 1) define the Hamiltonian or the potential respectively (--> function: potential) 2) minimize the potential ( I am using scipy.optimize.fmin to calculate the minimum of the potential) (2b) Expand the Hamiltonian around the minimum position. This step is not done in the code because it is only necessary for the purpose that one understand the physical background and why one have to do step 3) 3) Calculate the Hessian matrix of the potential, that means calculate the second derivatives of the potential at the point of the minimum position (--> numdifftools.Hessian) 4) Calculate the eigenvalues of the Hessian (-->scipy.linalg.eigvals) The problem can be solved analytically except of the calculation of the minima in step 2: Now I have the following problem: The Hessian matrix evaluated at the minimum position is wrong. What I checked so far: 1) The potential seems to be calculated correctly and the minimum position for two ions (N=2) is fine. 2) The Hessian matrix calculation works fine for several other functions. 3) The Hessian matrix is correct when I set the Coulomb interaction to zero (C=0) 4) The Hessian matrix is correct when I set C to some simple arbitrary integer numbers like 1, 5 ,8 .... 5) The Hesse matrix gives the correct number of the first part of the diagonal elements even with C=2.3e-28 ! But in that case the offdiagonal elements and the part of the diagonal element which refers to the second derivative in the Coulomb term is wrong! The offdiagonal term should be C/(abs(x_1-x_2))**3 = 2.6e-12 but it is 3.4e-24! 5) In principal Python can divide those numbers (2*2.3e-28)/(5.6e-6)**3 6) I played around with the accurateness of the calculation of the Hessian by varying the arguments numTerms, method and metOrder but that didn't change anything in my case. My source code is attached. Please keep in mind that I just started programing and Python is my first Language. I tried to do my best with the comments to make the code readable. Here it comes: # import a tool to use / as a symbol for normal division from __future__ import division #import system data import sys, os #Loading the required packages import scipy as sp import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import numdifftools as nd # and subpackages from scipy import * from scipy import linalg, optimize, constants #----------------------------------------------------------------------------------------- # Hamiltonian H=sum_i(p_i^2/(2m)+ 1/2 * m * w^2 x_i^2)+ sum_(i!=j)(a/|x_i-x_j|) #----------------------------------------------------------------------------------------- class classicalHamiltonian: def __init__(self): self.N = 2 #N is a scalar, it's the number of ions in the chain f = 1000000 #f is a scalar, it's the trap frequency self.w = 2*pi*f #w is a scalar, it's the angular velocity corresponding to the trap frequency self.C = (4*pi*constants.epsilon_0)**(-1)*constants.e**2 #C is a scalar, it's the Coulomb constant times the electronic charge in SI self.m = 39.96*1.66e-27 #m is the mass of a single trapped ion in the chain def potential(self, positionvector): #Defines the potential that is going to be minimized x= positionvector #x is an 1-d array (vector) of lenght N that contains the positions of the N ions w= self.w C= self.C m= self.m #First we consider the potential of the harmonic osszilator Vx = 0.5 * m * (w**2) * sum(x**2) for i in range(len(x)): for j in range(i+1, len(x)): Vx += C * (abs(x[i]-x[j]))**(-1) #then we add the coulomb interaction return Vx def initialposition(self): #Defines the initial position as an estimate for the minimize process N= self.N x_0 = r_[-(N-1)/2:(N-1)/2:N*1j] return x_0 def normal_modes(self, eigenvalues): #the computed eigenvalues of the matrix Vx are of the form (normal_modes)^2*m. m = self.m normal_modes = sqrt(eigenvalues/m) return normal_modes #C=(4*pi*constants.epsilon_0)**(-1)*constants.e**2 c=classicalHamiltonian() #print c.potential(array([-0.5, 0.5])) xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10) hessian = nd.Hessian(c.potential) eigenvalues = linalg.eigvals(hessian(xopt)) normal_modes = c.normal_modes(eigenvalues) I appreciate any kind of help. Thanks a lot in advance Nicolai ps: If you have any questions or need some more information please let me know.

3 2

loading and filtering data for matplotlib
by Resmi 26 Oct '10

26 Oct '10

Hi I'm new to python. I'm trying to plot a data set in matplotlib. The data is a mixture of strings and floats, has four columns (col[0],col[1],col[2],col[3]) and it looks like: R -2.29350 0.50340 0.480E-01 R -2.25903 0.50740 0.480E-01 SU -2.19457 0.16200 0.800E-01 SU -2.13237 0.14600 0.800E-01 What I want to do is to plot the SU filtered entries alone as col[1]:col[2] I wrote a script to read the data, how can I make x1 and SUx2 as arrays? Right now matplotlib only plots the last value obviously. Or is there a different/better way to slice just the SU entries and make them as a separate array? My actual data file has many more entries with more filters and I want to make a plot for each of the filters with col[1]:col[2] import pylab import numpy #Loading and reading data data=numpy.loadtxt('file.dat',dtype=[('filt','S4'),('x1','f8'),('x2','f8'),('x3','f8')]) filt=data['filt'] for i in range(0,len(filt),1): if filt[i] == 'SU': x1=data['x1'][i] SUx2=data['x2'][i] #print x1,SUx2 pylab.plot(x1,SUx2,'bs') pylab.show()

2 2

portable doctests despite floating points numbers
by Sébastien Barthélemy 26 Oct '10

26 Oct '10

Hello all, I use doctest for examples and tests in a program which relies heavily on numpy. As floating point calculations differs slightly across computers (32/64 bits), I have troubles writing portable doctests. The doctest documentation [1] advises to use numbers in the form int/2**n. This is quite restrictive. Using numpy.set_printoptions(suppress=True) also helps, but I still have problems with numbers around 0. On some platform a result prints as 0., on others as -0. Is there a workaround? [1] http://docs.python.org/library/doctest.html#warnings

2 3

speed of array vs matrix
by Citi, Luca 25 Oct '10

25 Oct '10

Hello, I have noticed a significant speed difference between the array and the matrix implementation of the dot product, especially for not-so-big matrices. For example: In [1]: import numpy as np In [2]: b = np.random.rand(104,1) In [3]: bm = np.mat(b) In [4]: a = np.random.rand(8, 104) In [5]: am = np.mat(a) In [6]: %timeit np.dot(a, b) 1000000 loops, best of 3: 1.74 us per loop In [7]: %timeit am * bm 100000 loops, best of 3: 6.38 us per loop The results for two different PCs (PC1 with windows/EPD6.2-2 and PC2 with ubuntu/numpy-1.3.0) and two different sizes are below: array matrix 8x104 * 104x1 PC1 1.74us 6.38us PC2 1.23us 5.85us 8x10 * 10x5 PC1 2.38us 7.55us PC2 1.56us 6.01us For bigger matrices the timings seem to asymptotically approach. Is it something worth trying to fix or should I just accept this as a fact and, when working with small matrices, stick to array? Thanks, Luca

3 3