Mailman 3 April 2008 - NumPy-Discussion

dot/tensordot limitations
by Anne Archibald 30 Apr '08

30 Apr '08

Timothy Hochberg has proposed a generalization of the matrix mechanism to support manipulating arrays of linear algebra objects. For example, one might have an array of matrices one wants to apply to an array of vectors, to yield an array of vectors: In [88]: A = np.repeat(np.eye(3)[np.newaxis,...],2,axis=0) In [89]: A Out[89]: array([[[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]], [[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]]]) In [90]: V = np.array([[1,0,0],[0,1,0]]) Currently, it is very clumsy to handle this kind of situation even with arrays, keeping track of dimensions by hand. For example if one wants to multiply A by V "elementwise", one cannot simply use dot: In [92]: np.dot(A,V.T) Out[92]: array([[[ 1., 0.], [ 0., 1.], [ 0., 0.]], [[ 1., 0.], [ 0., 1.], [ 0., 0.]]]) tensordot() suffers from the same problem. It arises because when you combine two multiindexed objects there are a number of ways you can do it: A: Treat all indices as different and form all pairwise products (outer product): In [93]: np.multiply.outer(A,V).shape Out[93]: (2, 3, 3, 2, 3) B: Contract the outer product on a pair of indices: In [98]: np.tensordot(A,V,axes=(-1,-1)).shape Out[98]: (2, 3, 2) C: Take the diagonal along a pair of indices: In [111]: np.diagonal(np.multiply.outer(A,V),axis1=0,axis2=3).shape Out[111]: (3, 3, 3, 2) What we want in this case is a combination of B and C: In [106]: np.diagonal(np.tensordot(A,V,axes=(-1,-1)),axis1=0,axis2=2).T.shape Out[106]: (2, 3) but it cannot be done without constructing a larger array and pulling out its diagonal. If this seems like an exotic thing to want to do, suppose instead we have two arrays of vectors, and we want to evaluate the array of dot products. None of no.dot, np.tensordot, or np.inner produce the results you want. You have to multiply elementwise and sum. (This also precludes automatic use of BLAS, if indeed tensordot uses BLAS.) Does it make sense to implement a generalized tensordot that can do this, or is it the sort of thing one should code by hand every time? Is there any way we can make it easier to do simple common operations like take the elementwise matrix-vector product of A and V? The more general issue, of making linear algebra natural by keeping track of which indices are for elementwise operations, and which should be used for dots (and how), is a whole other kettle of fish. I think for that someone should think hard about writing a full-featured multilinear algebra package (might as well keep track of coordinate transformations too while one was at it) if this is intended. Anne

4 5

Starting to work on runtime plugin system for plugin (automatic sse optimization, etc...)
by David Cournapeau 30 Apr '08

30 Apr '08

Hi, I've just started working on a prototype for a plugin system for numpy. The plugin aims at providing a framework for the following user cases: - runtime selection of blas/lapack/etc...: instead of harcoding in the binary one blas/lapack implementation, numpy could choose the SSE optimized if the CPU supports SSE, etc... - this could also be used for core numpy, for example ufuncs: if we want to start implementing some tight loop with aggressively optimized code (SSE, etc...), we could again ship with a default pure C implementation, and choose the best one at runtime. - we could even have a system to choose a different implementation (for example, right now, scipy is shipped with a slow fft for licensing issues mainly, and people installing fftw could then tell scipy to use fftw instead of the included one). Right now, the prototype does not do much, and only works for linux; I mainly focused on automatic generation of the plugin from a list of functions, and transparent use from numpy point of view. It provides the plugin api through pure function pointers, without the need for the user to be aware of it. For example, if you have an api with the following functions: void foo1(); int foo2(); int foo3(int); int foo4(double* , double*); int foo5(double* , double*, int); The current implementation would build the boilerplate to load those functions, etc... and you would just use those functions in numpy like the following: init_foo(); /* all functions are prefixed with npyw, for numpy wrapper */ npyw_foo1(); npyw_foo2(n); etc... The code can be found there: https://code.launchpad.net/~david-ar/+junk/numplug And some thinking (pretty low content for now): http://www.scipy.org/RuntimeOptimization cheers, David

6 16

Can anyone look at #759
by David Cournapeau 30 Apr '08

30 Apr '08

Hi, I have a patch for an issue which had bugged me for some time, but since I am not familiar with that part of numpy code, I would prefer someone checking I am not doing anything stupid before checking in, http://www.scipy.org/scipy/numpy/ticket/759 cheers, David

2 2

accuracy issues with numpy arrays?
by eli bressert 29 Apr '08

29 Apr '08

Hi, I'm writing a quick script to import a fits (astronomy) image that has very low values for each pixel. Mostly on the order of 10^-9. I have written a python script that attempts to take low values and put them in integer format. I basically do this by taking the mean of the 1000 lowest pixel values, excluding zeros, and dividing the rest of the image by that mean. Unfortunately, when I try this in practice, *all* of the values in the image are being treated as zeros. But, if I use a scipy.ndimage function, I get proper values. For example, I take the pixel that I know has the highest value and do x = scipy.ndimage.maximum(image) print x 1.7400700016878545e-05 The script is below. Thanks for the help. Eli import pyfits as p import scipy as s import scipy.ndimage as nd import numpy as n def flux2int(name): d = p.getdata(name) x,y = n.shape(d) l = x*y arr1 = n.array(d.reshape(x*y,1)) temp = n.unique(arr1[0]) # This is where the bug starts. All values are treated as zeros. Hence only one value remains, zero. arr1 = arr1.sort() arr1 = n.array(arr1) arr1 = n.array(arr1[s.where(arr1 >= temp)]) val = n.mean(arr1[0:1000]) d = d*(1.0/val) d = d.round() p.writeto(name[0,]+'fixed.fits',d,h)

2 1

Linear Interpolation 2
by Andrea Gavana 29 Apr '08

29 Apr '08

Hi All, as I didn't get any answer, I thought I might repost with some more analysis the trouble I am having with interp2d. The attached script produces 3 different results depending on the value of the parameter numColumns, i.e.: 1) numColumns = 20 Traceback (most recent call last): File "C:\MyProjects\LinearInterpolation.py", line 27, in <module> function = interp2d(xx, yy, z, kind="linear") File "C:\Python25\Lib\site-packages\scipy\interpolate\interpolate.py", line 113, in __init__ self.tck = fitpack.bisplrep(self.x, self.y, self.z, kx=kx, ky=ky, s=0.) File "C:\Python25\Lib\site-packages\scipy\interpolate\fitpack.py", line 702, in bisplrep tx,ty,nxest,nyest,wrk,lwrk1,lwrk2) ValueError: Invalid inputs. 2) numColumns = 200 Traceback (most recent call last): File "C:\MyProjects\LinearInterpolation.py", line 27, in <module> function = interp2d(xx, yy, z, kind="linear") File "C:\Python25\Lib\site-packages\scipy\interpolate\interpolate.py", line 113, in __init__ self.tck = fitpack.bisplrep(self.x, self.y, self.z, kx=kx, ky=ky, s=0.) File "C:\Python25\Lib\site-packages\scipy\interpolate\fitpack.py", line 702, in bisplrep tx,ty,nxest,nyest,wrk,lwrk1,lwrk2) MemoryError (MemoryError? With a 372x200 matrix???) 3) numColumns = 2000 Traceback (most recent call last): File "C:\MyProjects\LinearInterpolation.py", line 27, in <module> function = interp2d(xx, yy, z, kind="linear") File "C:\Python25\Lib\site-packages\scipy\interpolate\interpolate.py", line 113, in __init__ self.tck = fitpack.bisplrep(self.x, self.y, self.z, kx=kx, ky=ky, s=0.) File "C:\Python25\Lib\site-packages\scipy\interpolate\fitpack.py", line 702, in bisplrep tx,ty,nxest,nyest,wrk,lwrk1,lwrk2) OverflowError: long int too large to convert to int Is there anyone who could tell me, please, what I am doing wrong? This is on Windows XP 1GB RAM, 3.7 GHz, Python 2.5, scipy 0.6.0, numpy 1.0.4. Thank you for your suggestions. Andrea. "Imagination Is The Only Weapon In The War Against Reality." http://xoomer.alice.it/infinity77/

1 0

Re: [Numpy-discussion] Matrix Class [was numpy release]
by Timothy Hochberg 29 Apr '08

29 Apr '08

Chris.Barker wrote: > Alan G Isaac wrote: > > the cost of complexity should be justified by a gain in functionality. > > I don't think functionality is the right word here. the Matrix class(es) > is all about clean, convenient API, i.e. style, not functionality -- we > have all the functionality already, indeed we have it with plain old > arrays, so I think that's really beside the point. Not entirely, there's no good way do deal with arrays of matrices at present. This could be fixed by tweaking dot, but it could also be part of a reform of the matrix class. [CHOP] > > Timothy Hochberg wrote: > > 1. The matrices and arrays should become more alike if possible > > I'm not sure I agree -- the more alike they are, the less point there is > to having them. That's the best possible outcome. If some solution can be reached that naturally supports enough matrix operations on array, without significantly complexifying array, to satisfy the matrix users then that would be great. Less stuff to maintain, less stuff to learn, etc, etc. With that in mind, what is minimum amount of stuff that matrix should support: 1. Indexing along a single axis should produce a row or column vector as appropriate. 2. '*' should mean matrix multiplication. 3. '**' should mean matrix exponentiation. I suspect that this is less crucial than 2, but no more difficult. There's some other stuff that's less critical IMO (.H and .I for example) but feel free to yell at me if you think I'm mistaken. There's some other properties that a fix should have as well, in my opinion and in some cases in others opinions as well. 1. A single index into an array should yield a 1D object just as indexing an array does. This does not have to inconsistent with #1 above; Chris Barker proposed one solution. I'm not sold on the details of that solution, but I approve of the direction that it points to. [In general A[i][j] should equal A[i,j]. I know that fancy indexing breaks this; that was a mistake IMO, but that's a discussion for another day]. 2. It should be possible to embed both vectors and matrices naturally into arrays. This would allow natural and efficient implementation of, for example, rotating a 1000 3D vectors. One could spell that R * V, where R is a 2x2 matrix and V is a 1000x3 array, where the last axis is understood to be a vector. 3. I'm pretty certain there's a third thing I wanted to mention but it escapes me. It'll resurface at the most inopportune time.... Let's play with Chris Barker's RowVector and ColVector proposal for a moment. Let's suppose that we have four four classes: Scalar, RowVector, ColVector and Matrix. They're all pretty much the same as today's array class except that: 1. They are displayed a little differently so that you can tell them apart. 2. __mul__ and __pow__ are treated differently Let's consider __mul__: when a RowVector multiplied with ColVector, the dot product is computed. If, for example, the arrays are both 2D, the they are treated as arrays of vectors and the dot product of each pair of vectors is computed in turn: I think broadcasting should work correctly, but I haven't thought that all the way through yet. Ignoring broadcasting for a moment, the rules are: 1. (n)D Scalar * (n)D Scalar => (n)D Scalar (elementwise) 2. (n)D RowVector * (n)D ColVector => (n-1)D Scalar (dot product) 3. (n+1)D Matrix * (n)D ColVector => (n)D ColVector (matrix-vector product) 4. (n)D Matrix * n(D) Matrix => (n)D Matrix (matrix) product Other combinations would be an error. In principal you could do dyadic products, but I suspect we'd be better off using a separate function for that since most of the times that would just indicate a mistake. Note that in this example Scalar is playing the role of the present day array, which I've assumed has magically vanished from the scene somehow. This looks like it gets of most the way towards where we want to be. There are some issues though. For example, all of the sudden you want to different transpose operators; one that transposes matrices and flips colVectors to RowVectors and leaves Scalars alone, and another that manipulates the rest of the array structure. There is probably other stuff like that too, it doesn't look insurmountable to me right now, however. Still, I'm not sold on the whole RowVector/ColVector/Matrix approach. I have a gut feeling that we'd be better off with a single array class that was somehow tagged to indicate it's type. Also, I keep hoping to come up with an elegant generalization to this that turns arrays into quasi-tensor objects where all the matrix behavior falls out naturally. Sadly, attempts in this direction keep collapsing under there own weight but I'm still thinking about it. However, I do think that the RowVector/ColVector/Matrix approach is a step in the right direction and is certainly worth discussing to see where it leads. -tim > > > should share more of the same code base. > > Don't they share almost all the same code already? My understanding is > that they are arrays with a couple operators overloaded -- how much more > code could they share? > > Nevertheless, I'm interested to see what Tim comes up with. > > > 2. This doesn't support the concept of arrays of matrices/vectors, > > which is one thing I've always wanted out of the matrix class. For > > example it would be nice to be able to spell: 'A' is a 1D array of > > matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall, > > You could do this now with a 1-d object array, with a bunch of matrices > in it. > > > matrix multiply them together, yielding a 1D array of matrices > > (3D) overall. > > But this, or course, would fail. However, if we expand this scenario to > A and B being 1-d arrays of other arrays that may not all be the same > size (and thus you couldn't just use a 3-d array), and be able to run > various ufuncs on them, that would be cool! > > Bill Spotz wrote: > >> The fact that there is another way to get a row vector: M[i] is a > >> bonus. > > > > Except that the behavior of M[i] is one of the driving issues of the > > conversation. > > well, yes, but a bit of a red herring, I think -- it's syntactic sugar. > The issue, I think, is that even if you do: > > r = M[:,i] > > you can't then index that with a single index -- it's a matrix, NOT a > 1-d array. Again, I proposed disallowing M[i] altogether, as it is > ambiguous, but that really is gratuitous consistency. > > >> A vector is a 1-d array with some functionality overloaded to make it > >> convenient to do linear algebra. > > > > Again, I would argue for Vectors inheriting from Matrix. I would make > > the case based on code reuse and elegance, although these might be > > overridden by some other concern. > > I think that's a bit of an implementation detail -- we need to define > behavior that we want, and decide how best to implement that. i.e.: > > RowVector[i] -> scalar > RowVector.A -> 1-d array > (same for column) > > How do we want it pretty-printed? > > We need to focus on column vectors here -- It's easy to say: a row > vector is a 1-d array -- that's easy and clean. It's column vectors that > are tricky -- indeed they are a pain with regular old arrays. > > What about > > np.Matrix(<container of RowVector>) > > np.Matrix(<container of ColumnVector>) > > I'd say you get a matrix from each of those, but they would be > transposes of each-other -- another good reason for the Vector classes. > though maybe those should be spelled: > > np.row_stack and np.column_stack > > > I have generally thought about this in the context of, say, a > > Krylov-space iterative method, and what that type of interface would > > lead to the most readable code. > > Can you whip up a small example, starting with the current implementation? > > Bruce Southey wrote: > > Hi, > > I would like to use the matrix functionality but I, like others, get by > > without it. > > What does it lack that keeps you from using it? That's the key question? > > > +1 to Tim's and Nadav's comments. As Tim said, there should be seamless > > integration between concepts of vectors and matrices - after all there > > really no distinction between them in linear algebra. > > This is all about being able to index a "vector" with a single index -- > that's it, I think. Matlab handles this by special casing matrices that > have a single dimension of one. It also has an easy way to spell a > literal for column vector, though I suppose: > > np.Matrix((1,2,3,4,5)).T > > isn't so bad. > > > -2 for having rowvector and columnvector - both numpy and I should know > > the orientation, if not, I deserve garbage or an error. > > But HOW can you know the orientation? a 1-d array has no orientation -- > which is why the current version always returns a matrix when indexing. > > > 0 for vector - I don't see the need for it as a unique entity but > > understand the simplicity of saying vector. > > I don't think we can have vector without ColumnVector, though I suppose > a "Vector" can be either a (n,1) or a (1,n) matrix that allows single > indexing. > > By the way, maybe a column vector could be handy for doing array > broadcasting, as an cleaner way to do: > > x = np.arange(10) > y = np.arange(20).reshape((-1,1)) > > z = FunctionOfXandY(x,y) > > > Final note: > > If no one has an linear algebra examples, then we're wasting out time > with even this cheap talk. > > -Chris > > > -- > Christopher Barker, Ph.D. > Oceanographer > > Emergency Response Division > NOAA/NOS/OR&R (206) 526-6959 voice > 7600 Sand Point Way NE (206) 526-6329 fax > Seattle, WA 98115 (206) 526-6317 main reception > > Chris.Barker(a)noaa.gov > _______________________________________________ > Numpy-discussion mailing list > Numpy-discussion(a)scipy.org > http://projects.scipy.org/mailman/listinfo/numpy-discussion > -- . __ . |-\ . . tim.hochberg(a)ieee.org

5 9

Indexing a matrix with a scalar and ticket #707
by Charles R Harris 29 Apr '08

29 Apr '08

Yes, indeed. Ticket #707: numpy.array fails if the input is a list of matrixes (with more then one column). The subroutine discover_dimensions in arrayobject.c indexes a matrix with a scalar. It is a recursive routine and expects to find the next lower dimension as it recurses down into the matrix. This works fine with matrices with one column because the matrix row also has first dimension 1. Things don't do so well if the matrices have more than one column. I expect this is just the tip of a big pile of, um, stuff. The easy fix is to make indexing by scalars the same for matrices and arrays. So I think the best thing to do is provide that and figure out how to get the subarrays some other way. Chuck

1 2

leaving matrix indexing discussion
by Alan G Isaac 28 Apr '08

28 Apr '08

At this point, I am starting to feel that I have abused my interlocutors. So I will withdraw from the conversation. Apologies to those who feel I should have done so earlier. The core issues are, I think, on the discussion page. <URL:http://www.scipy.org/MatrixIndexing> I hope people will read it before presenting arguments. And improve it too. Cheers, Alan Isaac

3 2

setting rows of a numpy array
by wilson 28 Apr '08

28 Apr '08

hi all, i have a numpy array of floats whose rows i need to set. i am setting the each row using pixel values of each image in some folder.I wrote a class MyImage that has a field pixelarray which stores the pixels of that image as a numpy array.I made several MyImage instances and stored them in a list called myimagefileslist. i can get the pixels of each image as a numpy array using myimagefileslist[i]._pixelarray where _pixelarray is a field of MyImage class . now i defined a function that sets a row of a given array as below def putrow(myarray,inrow ,rownum): myarray[rownum,:]=inrow in my code i initially make an empty array myimagesdata=zeros((numofimgs,numofpixels)) then i set the rows using for i in range(numofimgs): pixarray=asfarray(myimagefileslist[i]._pixelarray) putrow(myimagesdata,pixarray,i) this gets me the myimagesdata array updated with the pixeldata of images from myimagefileslist.But i am wondering if there is a better way to do this.If anyone can suggest improvements it wd be helpful W

2 1

Distribution Functions Change Behavior
by Rich Shepard 28 Apr '08

28 Apr '08

AFter the extensive input from folks here last week, and my examination of alternatives, I accepted the visual appearance of Gaussian curves in our model. As I check the plots for data errors I find a behavior change when the x-axis length is 14 rather than 100, and I do not understand why. I would appreciate the ideas of the mathematically adroit here to help resolve the problem. However, if this is not the appropriate forum for such a discussion, please point me to the more appropriate one. Most of the plots cover a range of independent values of 0-100. Some are greater (e.g., 0-1,000 or 0-10,000). However, for pH the range is 0-14. The plots for pH do not consist of sigmoid curves (a 'Z-curve' from 0-7 and an 'S-curve' from 7-14); they are straight lines. The Gaussian curve does not approach zero on either end (0 or 14). I can provide the functions for these curves. I need to understand why the results change when the range of the independent variable is short and learn how to produce correct results. Rich -- Richard B. Shepard, Ph.D. | Integrity Credibility Applied Ecosystem Services, Inc. | Innovation <http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

4 5