Mailman 3 February 2008 - NumPy-Discussion

Faster array version of ndindex
by Jonathan Taylor 22 Feb '08

22 Feb '08

I was needing an array representation of ndindex since ndindex only gives an iterator but array(list(ndindex)) takes too long. There is prob some obvious way to do this I am missing but if not feel free to include this code which is much faster. In [252]: time a=np.array(list(np.ndindex(10,10,10,10,10,10))) CPU times: user 11.61 s, sys: 0.09 s, total: 11.70 s Wall time: 11.82 In [253]: time a=ndtuples(10,10,10,10,10,10) CPU times: user 0.32 s, sys: 0.21 s, total: 0.53 s Wall time: 0.60 def ndtuples(*dims): """Fast implementation of array(list(ndindex(*dims))).""" # Need a list because we will go through it in reverse popping # off the size of the last dimension. dims = list(dims) # N will keep track of the current length of the indices. N = dims.pop() # At the beginning the current list of indices just ranges over the # last dimension. cur = np.arange(N) cur = cur[:,np.newaxis] while dims != []: d = dims.pop() # This repeats the current set of indices d times. # e.g. [0,1,2] -> [0,1,2,0,1,2,...,0,1,2] cur = np.kron(np.ones((d,1)),cur) # This ranges over the new dimension and 'stretches' it by N. # e.g. [0,1,2] -> [0,0,...,0,1,1,...,1,2,2,...,2] front = np.arange(d).repeat(N)[:,np.newaxis] # This puts these two together. cur = np.column_stack((front,cur)) N *= d return cur

4 4

user-defined types mixed arithmetic
by Neal Becker 22 Feb '08

22 Feb '08

I have 2 user-defined types, and simple arithmetic working for the cases I registered with PyUFunc_RegisterLoopForType. I'd like to use automatic conversion to do mixed arithmetic between these 2 types. I did PyArray_RegisterCastFunc, and it seems this allows explicit conversion: >>> a array([(0,0), (1,0), (2,0), (3,0), (4,0), (5,0), (6,0), (7,0), (8,0), (9,0)], dtype=cmplx_int32) >>> b array([(0,0), (1,0), (2,0), (3,0), (4,0), (5,0), (6,0), (7,0), (8,0), (9,0)], dtype=cmplx_int64) >>> array (a,dtype=b.dtype) array([(0,0), (1,0), (2,0), (3,0), (4,0), (5,0), (6,0), (7,0), (8,0), (9,0)], dtype=cmplx_int64) >>> But mixed mode arithmetic gives: a+b TypeError: function not supported for these types, and can't coerce safely to supported types I've been trying to understand this 'coerce' without much luck. Any clues what I need to do here?

1 0

cmp_arg_types bug?
by Neal Becker 22 Feb '08

22 Feb '08

cmp_arg_types(int *arg1, int *arg2, int n) { while (n--) { if (PyArray_EquivTypenums(*arg1, *arg2)) continue; if (PyArray_CanCastSafely(*arg1, *arg2)) return -1; return 1; } return 0; } IIUC, if can cast (arg1, arg2), we never compare other args, just return -1. Shouldn't this be: cmp_arg_types(int *arg1, int *arg2, int n) { while (n--) { if (PyArray_EquivTypenums(*arg1, *arg2) or (PyArray_CanCastSafely(*arg1, *arg2)) continue; return 1; } return 0; }

1 1

ufunc for user-defined type
by Neal Becker 22 Feb '08

22 Feb '08

Now that I have my user-defined type, I want to add some funcs. According to the numpy book, I need to use: PyUFunc_RegisterLoopForType The book says I first need a ufunc. The only way I see to create one is PyUFunc_FromFuncAndData. Is the the correct procedure? I wonder, because PyUFunc_FromFuncAndData requires 'type', and 'type' is char, but user-defined types start at 256, which doesn't fit in a char, which gives a compile warning.

2 2

change memmap.sync function
by Christopher Burns 22 Feb '08

22 Feb '08

Would anyone oppose deprecating the memmap.sync function and replacing it with memmap.flush? This would match python's mmap module, and I think be more intuitive. -- Christopher Burns, Software Engineer Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

2 1

finding eigenvectors etc
by devnew＠gmail.com 21 Feb '08

21 Feb '08

hi i was calculating eigenvalues and eigenvectors for a covariancematrix using numpy adjfaces=matrix(adjarr) faces_trans=adjfaces.transpose() covarmat=adjfaces*faces_trans evalues,evect=eigh(covarmat) for a sample covarmat like [[ 1.69365981e+13 , -5.44960784e+12, -9.00346400e+12 , -2.48352625e +12] [ -5.44960784e+12, 5.08860660e+12, -8.67539205e+11 , 1.22854045e +12] [ -9.00346400e+12, -8.67539205e+11, 1.78184943e+13 ,-7.94749110e +12] [ -2.48352625e+12 , 1.22854045e+12, -7.94749110e+12 , 9.20247690e +12]] i get these evalues [ 3.84433376e-03, 4.17099934e+12 , 1.71771364e+13 , 2.76980401e+13] evect [[ 0.5 -0.04330262 0.60041892 -0.62259297] [ 0.5 -0.78034307 -0.35933516 0.10928372] [ 0.5 0.25371931 0.3700265 0.74074753] [ 0.5 0.56992638 -0.61111026 -0.22743827]] what bothers me is that for the same covarmat i get a different set of eigenvectors and eigenvalues when i use java library Jama's methods Matrix faceM = new Matrix(faces, nrfaces,length); Matrix faceM_transpose = faceM.transpose(); Matrix covarM = faceM.times(faceM_transpose); EigenvalueDecomposition E = covarM.eig(); double[] eigValue = diag(E.getD().getArray()); double[][] eigVector = E.getV().getArray(); here the eigValue= [-6.835301757686207E-4, 4.170999335736721E12, 1.7177136443134865E13, 2.7698040117669414E13] and eigVector [ [0.5, -0.04330262221379265, 0.6004189175979487, 0.6225929700052174], [0.5, -0.7803430730840767, -0.3593351608695496, -0.10928371540423852], [0.49999999999999994, 0.2537193127299541, 0.370026504572483, -0.7407475253159538], [0.49999999999999994, 0.5699263825679145, -0.6111102613008821, 0.22743827071497524] ] I am quite confused bythis difference in results ..the first element in eigValue is different and also the signs in last column of eigVectors are diff..can someone tell me why this happens? thanks dn

6 10

Rename record array fields
by Sameer DCosta 20 Feb '08

20 Feb '08

Is there a way to rename record array fields without making a copy of the whole record array? Thanks in advance for your replies. Sameer

3 2

bug in numpy.histogram?
by John Hunter 20 Feb '08

20 Feb '08

We recently deprecated matplotlib.mlab.hist, and I am now hitting a bug in numpy's historgram, which appears to be caused by the use of "any" that does not exist in the namespace. Small patch attached. The example below exposes the bug: Python 2.4.2 (#1, Feb 23 2006, 12:48:31) Type "copyright", "credits" or "license" for more information. IPython 0.8.3.svn.r2876 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [1]: import numpy as np In [2]: np.__file__ Out[2]: '/home/titan/johnh/dev/lib/python2.4/site-packages/numpy/__init__.pyc' In [3]: np.__version__ Out[3]: '1.0.5.dev4814' In [4]: x = np.random.randn(100) In [5]: bins = np.linspace(x.min(), x.max(), 40) In [6]: y = np.histogram(x, bins=bins) ------------------------------------------------------------ Traceback (most recent call last): File "<ipython console>", line 1, in ? File "/home/titan/johnh/dev/lib/python2.4/site-packages/numpy/lib/function_base.py", line 155, in histogram if(any(bins[1:]-bins[:-1] < 0)): NameError: global name 'any' is not defined

3 3

numpy record array segfault
by Sameer DCosta 20 Feb '08

20 Feb '08

Hi, I'm getting a segfault when using python objects with record arrays. The code (below) basically assigns a single datetime object to a slice of a column in the record array and then python segfaults as soon as I try to access those array values. I'm using the latest svn version of numpy compiled with gcc 3.4.1 on Solaris (intel). Any ideas why this is happening? A reference counting problem maybe? Thanks for taking a look. Sameer johnh@flag:~> gcc --version gcc (GCC) 3.4.1 Copyright (C) 2004 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. johnh@flag:~> uname -a SunOS flag 5.10 Generic_118855-15 i86pc i386 i86pc johnh@flag:~> python -V Python 2.4.2 johnh@flag:~> cat tmp.py import datetime import numpy as np print np.__version__ def myfunc(N): newrec = np.empty(N, dtype=[('date', '|O4'), ('age', int)]) newrec['date'] = datetime.date(2002,1,1) newrec['age'] = 22 newrec['date'][1:12] = datetime.date(2003,1,1) return newrec.view(np.recarray) if __name__=='__main__': newrec = myfunc(29) print newrec['date'] johnh@flag:~> python ~/tmp.py 1.0.5.dev4812 Segmentation Fault (core dumped)

2 2

mixed mode arithmetic
by Neal Becker 20 Feb '08

20 Feb '08

I've been browsing the numpy source. I'm wondering about mixed-mode arithmetic on arrays. I believe the way numpy handles this is that it never does mixed arithmetic, but instead converts arrays to a common type. Arguably, that might be efficient for a mix of say, double and float. Maybe not. But for a mix of complex and a scalar type (say, CDouble * Double), it's clearly suboptimal in efficiency. So, do I understand this correctly? If so, is that something we should improve?

1 0

NumPy-Discussion February 2008