Mailman 3 February 2008 - NumPy-Discussion

image to array doubt
by devnew＠gmail.com 29 Feb '08

29 Feb '08

hi i came across a codebase by rice univ people..in that there are some functions for conversion btw image and vectors 1.<code> def image_to_vector(self, filename): try: im = Image.open(filename) except IOError: print 'couldn\'t load ' + filename sys.exit(1) self.im_size = im.size a = numpy.array(im.getdata()) return a - 128 </code> what is the purpose of a-128 here? i couldn't quite understand 2.<code> def vector_to_image(self, v, filename): v.shape = (-1,) a, b = min(v), max(v) span = max(abs(b), abs(a)) im = Image.new('L', self.im_size) im.putdata((v * 127. / span) + 128) im.save(filename) </code> and why the calculations of a,b span? why do you need the extra maths ops inside im.putdata() I couldn't quite follow it can someone explain this? thanks D

3 3

ANN: sympycore version 0.1 released
by Pearu Peterson 29 Feb '08

29 Feb '08

We are proud to present a new Python package: sympycore - an efficient pure Python Computer Algebra System Sympycore is available for download from http://sympycore.googlecode.com/ Sympycore is released under the New BSD License. Sympycore provides efficient data structures for representing symbolic expressions and methods to manipulate them. Sympycore uses a very clear algebra oriented design that can be easily extended. Sympycore is a pure Python package with no external dependencies, it requires Python version 2.5 or higher to run. Sympycore uses Mpmath for fast arbitrary-precision floating-point arithmetic that is included into sympycore package. Sympycore is to our knowledge the most efficient pure Python implementation of a Computer Algebra System. Its speed is comparable to Computer Algebra Systems implemented in compiled languages. Some comparison benchmarks are available in * http://code.google.com/p/sympycore/wiki/Performance * http://code.google.com/p/sympycore/wiki/PerformanceHistory and it is our aim to continue seeking for more efficient ways to manipulate symbolic expressions: http://cens.ioc.ee/~pearu/sympycore_bench/ Sympycore version 0.1 provides the following features: * symbolic arithmetic operations * basic expression manipulation methods: expanding, substituting, and pattern matching. * primitive algebra to represent unevaluated symbolic expressions * calculus algebra of symbolic expressions, unevaluated elementary functions, differentiation and polynomial integration methods * univariate and multivariate polynomial rings * matrix rings * expressions with physical units * SympyCore User's Guide and API Docs are available online. Take a look at the demo for sympycore 0.1 release: http://sympycore.googlecode.com/svn/trunk/doc/html/demo0_1.html However, one should be aware that sympycore does not implement many features that other Computer Algebra Systems do. The version number 0.1 speaks for itself:) Sympycore is inspired by many attempts to implement CAS for Python and it is created to fix SymPy performance and robustness issues. Sympycore does not yet have nearly as many features as SymPy. Our goal is to work on in direction of merging the efforts with the SymPy project in the near future. Enjoy! * Pearu Peterson * Fredrik Johansson Acknowledgments: * The work of Pearu Peterson on the SympyCore project is supported by a Center of Excellence grant from the Norwegian Research Council to Center for Biomedical Computing at Simula Research Laboratory.

1 0

I,J,K Coordinates from Cell ID
by Andrea Gavana 29 Feb '08

29 Feb '08

Hi All, I have some problems in figuring out a solution for an issue I am trying to solve. I have a 3D grid of dimension Nx, Ny, Nz; for every cell of this grid, I calculate the cell centroids (with the cell coordinates x, y, and z) and then I try to find which cell centroid is the closest to a specified point in 3D (which I supply). I still haven't figured out how to do this, even if I have some ideas (and suggestions for this problem are welcome :-D ). But the problem is another one. When I find the closest centroid, I know only the cell ID of this centroid. The cell ID is (usually) defined as: ID = (K-1)*Nx*Ny + (J-1)*Nx + I - 1 Where I, J, K are the cell indexes. Now, the problem is, how can I calculate back the I, J, K indexes knowing only the cell ID? I am trying to solve this using numpy (as my grid is stored using a numpy array), but it's something akin to the Matlab function ind2sub... Thank you for your suggestions. Andrea. "Imagination Is The Only Weapon In The War Against Reality." http://xoomer.alice.it/infinity77/

2 2

arrays of matrices
by Geoffrey Irving 28 Feb '08

28 Feb '08

Hello, I have a large number of points (shape (n,3)), and a matching number of 3x3 matrices (shape (n,3,3)), and I want to compute the product of each matrix times the corresponding point. I can't see a way to do this operation with dot or tensordot, since these routines either sum across an index or treat it as independent between the two arguments. Is this case, I can use the fact that 3 is small to do it, but is there a clean way for handle this kind of "array of matrix" situation in general? Thanks, Geoffrey

2 4

Handling of numpy.power(0, <something>)
by Stuart Brorson 28 Feb '08

28 Feb '08

I have been poking at the limits of NumPy's handling of powers of zero. I find some results which are disturbing, at least to me. Here they are: In [67]: A = numpy.array([0, 0, 0]) In [68]: B = numpy.array([-1, 0, 1+1j]) In [69]: numpy.power(A, B) Out[69]: array([ 0.+0.j, 1.+0.j, 0.+0.j]) IMO, the answers should be [Inf, NaN, and NaN]. The reasons: ** 0^-1 is 1/0, which is infinity. Not much argument here, I would think. ** 0^0: This is problematic. People smarter than I have argued for both NaN and for 1, although I understand that 1 is the preferred value nowadays. If the NumPy gurus also think so, then I buy it. ** 0^(x+y*i): This one is tricky; please bear with me and I'll walk through the reason it should be NaN. In general, one can write a^(x+y*i) = (r exp(i*theta))^(x+y*i) where r, theta, x, and y are all reals. Then, this expression can be rearranged as: (r^x) * (r^i*y) * exp(i*theta*(x+y*i)) = (r^x) * (r^i*y) * exp(i*theta*x) * exp(-theta*y) Now consider what happens to each term if r = 0. -- r^x is either 0^<positive> = 1, or 0^<negative> = Inf. -- r^(i*y) = exp(i*y*ln(r)). If y != 0 (i.e. complex power), then taking the ln of r = 0 is -Inf. But what's exp(i*-Inf)? It's probably NaN, since nothing else makes sense. Note that if y == 0 (real power), then this term is still NaN (y*ln(r) = 0*ln(0) = Nan). However, by convention, 0^<real> is something other than NaN. -- exp(i*theta*x) is just a complex number. -- exp(-theta*y) is just a real number. Therefore, for 0^<complex> we have Inf * NaN * <complex> * <real>, which is NaN. Another observation to chew on. I know NumPy != Matlab, but FWIW, here's what Matlab says about these values: >> A = [0, 0, 0] A = 0 0 0 >> B = [-1, 0, 1+1*i] B = -1.0000 0 1.0000 + 1.0000i >> A .^ B ans = Inf 1.0000 NaN + NaNi Any reactions to this? Does NumPy just make library calls when computing power, or does it do any trapping of corner cases? And should the returns from power conform to the above suggestions? Regards, Stuart Brorson Interactive Supercomputing, inc. 135 Beaver Street | Waltham | MA | 02452 | USA http://www.interactivesupercomputing.com/

6 9

numpy.test() fails if it runs after scipy.test()
by Loïc BERTHE 28 Feb '08

28 Feb '08

Hi, I've got one failure if I run numpy.test() after running scipy.test() : ====================================================================== ERROR: Ticket #396 ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/loic/tmp/Linux/lib/python2.5/site-packages/numpy/core/tests/test_regression.py", line 600, in check_poly1d_nan_roots self.failUnlessRaises(N.linalg.LinAlgError,getattr,p,"r") File "/home/loic/tmp/Linux/lib/python2.5/unittest.py", line 320, in failUnlessRaises callableObj(*args, **kwargs) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/numpy/lib/polynomial.py", line 623, in __getattr__ return roots(self.coeffs) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/numpy/lib/polynomial.py", line 124, in roots roots = _eigvals(A) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/numpy/lib/polynomial.py", line 37, in _eigvals return eigvals(arg) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/scipy/linalg/decomp.py", line 378, in eigvals return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/scipy/linalg/decomp.py", line 128, in eig a1 = asarray_chkfinite(a) File "/home/loic/tmp/Linux/lib/python2.5/site-packages/numpy/lib/function_base.py", line 398, in asarray_chkfinite raise ValueError, "array must not contain infs or NaNs" ValueError: array must not contain infs or NaNs But I've got no error if I begin with numpy test. I've seen that this Ticket #396 seems closed in Trac, should I reopen it ? for more information, I've attached the results of scipy.test and numpy.test . Regards, -- LB

3 2

loadtxt broken if file does not end in newline
by Lisandro Dalcin 28 Feb '08

28 Feb '08

Dear all, I believe the current 'loadtxt' function is broken if file does not end in newline. The problem is at the last line of this fragment: for i,line in enumerate(fh): if i<skiprows: continue line = line[:line.find(comments)].strip() This is because `str.find()` (yet again!) returns -1 (minus one) on failure. Therefore, if a input ASCII file is missing the newline at the very end, the last character is chopped. Can someone file a bug or perhaps fix it? I'm very, very busy as to do any of that right now, sorry. Regards, -- Lisandro Dalcín --------------- Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC) Instituto de Desarrollo Tecnológico para la Industria Química (INTEC) Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) PTLC - Güemes 3450, (3000) Santa Fe, Argentina Tel/Fax: +54-(0)342-451.1594

7 17

ANN: Enthought Python Distribution - Beta
by Travis Vaught 28 Feb '08

28 Feb '08

Greetings, Enthought is very excited about our pending wide-release of the Enthought Python Distribution (EPD). After much effort, we finally think we're close to the first non-beta release. As one more quality check, we'd love to impose on you guys one more time to try out a just- minted beta release for Windows (EPD 2.5.2001_beta1) and give us some feedback. Any major problems will, of course, be fixed for the next release, but we're open to any suggestions for improvement for future releases as well. http://www.enthought.com/epd For those of you unfamiliar with EPD, it's a "kitchen-sink-included" distribution of Python with over 60 additional tools and libraries. It's bundled into a nice MSI installer on Windows and includes NumPy, SciPy, IPython, 2D and 3D visualization, database adapters and a lot of other tools right out of the box. We'll have support for RedHat and Mac OS X in a general release very soon. For academic, non-profit or hobbyist use, EPD is, and will remain, free. We are charging an annual subscription for commercial and governmental access to downloads and updates of EPD. Downloaded files may be used indefinitely past the subscription term. You are welcome to try out the beta indefinitely, regardless of your commercial/non- commercial persuasion. When the final (non-beta) version is released, commercial folks can try it for 30 days. You can check out the license terms (http://www.enthought.com/products/epdlicense.php) if you're interested in the details. EPD is compelling because it solves a lingering packaging and distribution problem, but also because of the libraries which it includes. We owe many folks on this list a debt of gratitude for their work on some really great tools. So, thanks ... and enjoy! Best Regards, Travis N. Vaught

4 3

Trouble With MaskedArray and Shared Masks
by Alexander Michael 27 Feb '08

27 Feb '08

I'm having trouble with MaskedArray's _sharedmask flag. I would like to create a sub-view of a MaskedArray, fill it, and have the modifications reflected in the original array. This works with regular ndarrays, but only works with MaskedArrays if _sharedmask is set to False. Here's an example: >>> a = numpy.ma.MaskedArray( ... data=numpy.zeros((4,5), dtype=float), ... mask=numpy.ones((4,5), dtype=numpy.ma.MaskType), ... fill_value=0.0 ... ) >>> sub_a = a[:2,:3] >>> sub_a[0,0] = 1.0 >>> print sub_a [[1.0 -- --] [-- -- --]] >>> print a [[-- -- -- -- --] [-- -- -- -- --] [-- -- -- -- --] [-- -- -- -- --]] >>> print a.data [[ 1. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]] The data array receives the new value, but the mask array does not. >>> a._sharedmask = False >>> sub_a = a[:2,:3] >>> sub_a[0,0] = 1.0 >>> print sub_a [[1.0 -- --] [-- -- --]] >>> print a [[1.0 -- -- -- --] [-- -- -- -- --] [-- -- -- -- --] [-- -- -- -- --]] This sort of (for me) unexpected behavior extends to other ways I've been using numpy arrays as well: a[:] = 1.0 (set to constant); a[:] = b (copy into); a[:5] = a[-5:] (rotating copy), etc. I wasn't seeing this behavior before because I was working on an array that had already been sliced and therefore "unshared", which caused a good deal of confusion for me when I started working on an array that wasn't the product of slicing. All of this leads me to some questions. What is the rational for initializing a new MaskedArray with _sharedmask=True when its mask isn't (actively) being shared yet? Is there a better way to say: "a=MaskedArray(...); a._sharedmask=False" that does not require touching a "private" attribute? Or am I going about this all wrong? What's the correct MaskedArray idioms for these actions that doesn't cause a new mask to be created? Thanks! Alex

2 3

Optimize speed of for loop using numpy
by Trond Kristiansen 27 Feb '08

27 Feb '08

Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this: tmp1=zeros((eta,xi),float) tmp2=zeros((eta,xi),float) tmp1=tmp1+10000 tmp2=tmp2+10000 for i in range(xi): for j in range(eta): for k in range(s): if z_r[k,j,i] < depth: if tmp1[j,i]==10000: if (depth - z_r[k,j,i]) <= (depth - z_r[0,j,i]): tmp1[j,i]=k else: if (depth - z_r[k,j,i]) <= (depth - z_r[int(tmp1[j,i]),j,i]): tmp1[j,i]=k elif z_r[k,j,i] >= depth: if tmp2[j,i]==10000: if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[s-1,j,i]) : tmp2[j,i]=k else: if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[int(tmp2[j,i]),j,i]) : tmp2[j,i]=k Not very impressive. My problem is that I can not find any numpy functions that actually can do the tests I do for each k value in the arrays. I need to identify the position (i,j) of k-values that meet the specified requirement. There are way too many if-else tests here as well. The scripts takes about 10 seconds to run (for 238MB input file), but these arrays are read from netCDF files and can be much larger and easily grow to enormous dimensions. It is therefore crucial for me to speed things up. Hope some of you can help. I really appreciate all feedback on this. I am just a rooky to numpy. Cheers and thanks for all help, Trond

6 8