Mailman 3 February 2009 - NumPy-Discussion

inplace matrix multiplication
by Frédéric Bastien April 29, 2009

April 29, 2009

Hi, I would like to know how I can make a call to the blas function gemm in numpy. I need a multiply and accumulate for matrix and I don't want to allocate a new matrix each time I do it. thanks for your time Frederic Bastien

4 6

Buildbot not building?
by Pauli Virtanen March 19, 2009

March 19, 2009

Hi, It seems that the buildbot.scipy.org is not picking up the changes in Numpy trunk. I'd guess this could be some issue with SVNPoller. At least it doesn't preserve states across buildmaster restarts, so replacing it with the following might help: {{{ import os from buildbot.changes.svnpoller import SVNPoller class PersistentSVNPoller(SVNPoller): persist_file = "svnpoll.rev" def __init__(self, *a, **kw): self.persist_file = kw.pop('persist_file', self.persist_file) SVNPoller.__init__(self, *a, **kw) if os.path.isfile(self.persist_file): f = open(self.persist_file, 'r') try: self.last_change = int(f.read()) except ValueError: pass finally: f.close() def get_new_logentries(self, *a, **kw): r = SVNPoller.get_new_logentries(self, *a, **kw) f = open(self.persist_file, 'w') try: f.write("%d" % self.last_change) finally: f.close() return r }}} -- Pauli Virtanen

2 8

Slicing/selection in multiple dimensions simultaneously
by Mike Ressler March 17, 2009

March 17, 2009

The following seems to be a wart: is it expected? Set up a 10x10 array and some indexing arrays: a=arange(100) a.shape=(10,10) q=array([0,2,4,6,8]) r=array([0,5]) Suppose I want to extract only the "even" numbered rows from a - then print a[q,:] <works - output deleted> Every fifth column: print a[:,r] <works - output deleted> Only the even rows of every fifth column: print a[q,r] --------------------------------------------------------------------------- <type 'exceptions.ValueError'> Traceback (most recent call last) /.../.../.../<ipython console> in <module>() <type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape But, this works: print a[q,:][:,r] [[ 0 5] [20 25] [40 45] [60 65] [80 85]] So why does the a[q,r] form have problems? Thanks for your insights. Mike -- mike.ressler(a)alum.mit.edu

8 16

inplace dot products
by Olivier Grisel March 16, 2009

March 16, 2009

Hi numpist people, I discovered the ufunc and there ability to compute the results on preallocated arrays: >>> a = arange(10, dtype=float32) >>> b = arange(10, dtype=float32) + 1 >>> c = add(a, b, a) >>> c is a True >>> a array([ 1., 3., 5., 7., 9., 11., 13., 15., 17., 19.], dtype=float32) My questions is : is there a way to have an equivalent for the dot product operation: I want atlas to build my dot products without allocating a temporary array and reuse a preallocated array of results. Suppose I have: >>> results = array((10, 3), dtype=float32) >>> W = arange(6, dtype=float32).reshape((2, 3)) >>> x = arange(20, dtype=float32).reshape((10, 2)) What I want is the equivalent of the following without the intermediate call to malloc: >>> results[:] = dot(x, W) Any idea? I tried to introspect the various docstring of numpy core modules but I could not get any lead. -- Olivier

4 9

image processing using numpy-scipy?
by Prashant Saxena March 12, 2009

March 12, 2009

Hi, This a little wiered problem. I am having a black and white image. (black background) Entire image is filled with noisy white patterns of different size and shape. I need to fill the white patches if there area is more then given one. Logically this could possible to use a quickfill algorithm for every pixel and if the filled area generated by that pixel is more then given area(in pixel) then paint that patch (filled area) with black. I have read the docs of PIL but there is no function for this. Can I use numpy-scipy for the matter? The image size is 1K. Any help would be greatly apperciated. Regards Prashant Connect with friends all over the world. Get Yahoo! India Messenger at http://in.messenger.yahoo.com/?wm=n/

4 5

Automatic differentiation (was Re: second-order gradient)
by Rob Clewley March 11, 2009

March 11, 2009

> Maybe we should focus on writing a decent 'deriv' function then. I > know Konrad Hinsen's Scientific had a derivatives package > (Scientific.Functions.Derivatives) that implemented automatic > differentiation: > > http://en.wikipedia.org/wiki/Automatic_differentiation That would be great, but wouldn't that be best suited as a utility requiring Sympy? You'll want to take advantage of all sorts of symbolic classes, especially for any source code transformation approach. IMO Hinsen's implementation isn't a very efficient or attractive solution to AD given the great existing C/C++ codes out there. Maybe we should be looking to provide a python interface to an existing open source package such as ADOL-C, but I'm all in favour of a new pure python approach too. What would be perfect is to have a single interface to a python AD package that would support a faster implementation if the user wished to install a C/C++ package, otherwise would default to a pure python equivalent. -Rob

4 4

N-D array interface page is out of date
by Andrew Straw March 9, 2009

March 9, 2009

Hi, I just noticed that the N-D array interface page is outdated and doesn't mention the buffer interface that is standard with Python 2.6 and Python 3.0: http://numpy.scipy.org/array_interface.shtml This page is linked to from http://numpy.scipy.org/ I suggest, at the minimum, modifying the page with really annoying blinking red letters at the top (or other suitable warning) that this is deprecated at that people should use http://www.python.org/dev/peps/pep-3118/ instead.

4 8

strange multiplication behavior with numpy.float64 and ndarray subclass
by Darren Dale March 8, 2009

March 8, 2009

I have a simple test script here that multiplies an ndarray subclass with another number. Can anyone help me understand why each of these combinations returns a new instance of MyArray: mine = MyArray() print type(np.float32(1)*mine) print type(mine*np.float32(1)) print type(mine*np.float64(1)) print type(1*mine) print type(mine*1) but this one returns a np.float64 instance? print type(np.float64(1)*mine) Here is the full script: import numpy as np class MyArray(np.ndarray): __array_priority__ = 20 def __new__(cls): return np.asarray(1).view(cls).copy() def __repr__(self): return 'my_array' __str__ = __repr__ def __mul__(self, other): return super(MyArray, self).__mul__(other) def __rmul__(self, other): return super(MyArray, self).__rmul__(other) mine = MyArray() print type(np.float32(1)*mine) print type(mine*np.float32(1)) print type(mine*np.float64(1)) print type(1*mine) print type(mine*1) print type(np.float64(1)*mine) Thanks, Darren

3 9

numpy.fix and subclasses
by Darren Dale March 8, 2009

March 8, 2009

I've been finding some numpy functions that could maybe be improved to work better with ndarray subclasses. For example: def fix(x, y=None): x = nx.asanyarray(x) if y is None: y = nx.zeros_like(x) y1 = nx.floor(x) y2 = nx.ceil(x) y[...] = nx.where(x >= 0, y1, y2) return y This implementation is a problematic for subclasses, since it does not allow metadata to propagate using the usual ufunc machinery of __array_wrap__, like ceil and floor do. nx.zeros_like does yield another instance of type(x), but y does not get x's metadata (such as units or a mask). Would it be possible to do something like: if y is None: y = x*0 "where" is another function that could maybe be improved to work with the rules established by array_priority, but I'm a lousy C programmer and I haven't actually looked into how this would work. If "where" respected array_priority, fix could be implemented as: def fix(x, y=None): x = nx.asanyarray(x) y1 = nx.floor(x) y2 = nx.ceil(x) if y is None: return nx.where(x >= 0, y1, y2) y[...] = nx.where(x >= 0, y1, y2) return y Darren

1 2

possible bug: __array_wrap__ is not called during arithmetic operations in some cases
by Darren Dale March 8, 2009

March 8, 2009

Does anyone know why __array_wrap__ is not called for subclasses during arithmetic operations where an iterable like a list or tuple appears to the right of the subclass? When I do "mine*[1,2,3]", array_wrap is not called and I get an ndarray instead of a MyArray. "[1,2,3]*mine" is fine, as is "mine*array([1,2,3])". I see the same issue with division, addition, etc. Here is a demonstration, observed with svn 6456: import numpy as np class MyArray(np.ndarray): __array_priority__ = 20 def __new__(cls): return np.asarray(1).view(cls).copy() def __array_wrap__(self, obj, context=None): print 'array wrap:', self, obj, context return obj.view(type(self)) def __str__(self): return 'MyArray(%s)'%super(MyArray,self).__str__() mine = MyArray() print 3*mine print mine*3 print [1,2,3]*mine print mine*[1,2,3] print print 3/mine print mine/3 print [1,2,3]*mine print mine*[1,2,3] print print 3+mine print mine+3 print [1,2,3]+mine print mine+[1,2,3] print print 3/mine print mine/3 print [1,2,3]/mine print mine/[1,2,3] print print 3**mine print mine**3 print [1,2,3]**mine print mine**[1,2,3]

5 15