Re: [Numpy-discussion] GPU Numpy

2009-08-06 Thread Sturla Molden
Olivier Grisel wrote:
 As usual, MS reinvents the wheel with DirectX Compute but vendors such
 as AMD and nvidia propose both the OpenCL API +runtime binaries for
 windows and their DirectX Compute counterpart, based on mostly the
 same underlying implementation, e.g. CUDA in nvidia's case.

   
Here is a DirectX Compute tutorial I found:

http://www.gamedev.net/community/forums/topic.asp?topic_id=516043

It pretty much says all we need to know. I am not investing any of my 
time learning that shitty API. Period.

Lets just hope OpenCL makes it to Windows without Microsoft breaking it 
for security reasons (as they did with OpenGL).



Sturla



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] maximum value and corresponding index

2009-08-06 Thread Hans Meine
On Wednesday 05 August 2009 22:06:03 David Goldsmith wrote:
 But you can cheat and put them on one line (if that's all you're after):
  x = np.array([1, 2, 3])
  maxi = x.argmax(); maxv = x[maxi]

Is there any reason not to put this as a convenience function into numpy?
It is needed so frequently, and it's a shame that the shortest solution 
traverses the array twice (the above is usually longer due to variable names, 
and/or 1 dimensions).

def give_me_a_good_name(array):
pos = array.argmax()
val = array[unravel_index(pos)]
return pos, val

The name should not be too long, maybe findmax?  I don't see good 
predecessors to learn from ATM.

OTOH, a minmax() would be more pressing, since there is no shortcut yet 
AFAICS.

Ciao,
  Hans
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] GPU Numpy

2009-08-06 Thread Romain Brette
Charles R Harris charlesr.harris at gmail.com writes:

 
 What sort of functionality are you looking for? It could be you could slip in
a small mod that would do what you want. In the larger picture, the use of GPUs
has been discussed on the list several times going back at least a year. The
main problems with using GPUs were that CUDA was only available for nvidia video
cards and there didn't seem to be any hope for a CUDA version of LAPACK.Chuck
 

So for our project what we need is:
* element-wise operations on vectors (arithmetical, exp/log, exponentiation)
* same but on views (x[2:7])
* assignment (x[:]=2*y)
* boolean operations on vectors (x2.5) and the nonzero() method
* possibly, multiplying a N*M matrix by an M*M matrix, where N is large and M is
small (but this could be done with vector operations).
* random number generation would be great too (gpurand(N))

What is very important to me is that the syntax be the same as with normal
arrays, so that you could easily switch the GPU on/off (depending on whether a
GPU was detected).

Cheers
Romain


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] GPU Numpy

2009-08-06 Thread Romain Brette
Ian Mallett geometrian at gmail.com writes:

 
 
 On Wed, Aug 5, 2009 at 11:34 AM, Charles R Harris charlesr.harris at
gmail.com wrote:
 
 
 
 It could be you could slip in a small mod that would do what you want.
 I'll help, if you want.  I'm good with GPUs, and I'd appreciate the numerical
power it would afford.

That would be great actually if we could gather a little team! Anyone else
interested?

As Trevor said, OpenCL could be a better choice than Cuda.

By the way, there is a Matlab toolbox that seems to have similar functionality:
http://www.accelereyes.com/

Cheers
Romain


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] GPU Numpy

2009-08-06 Thread Romain Brette
David Warde-Farley dwf at cs.toronto.edu writes:
 It did inspire some of our colleagues in Montreal to create this,  
 though:
 
   http://code.google.com/p/cuda-ndarray/
 
 I gather it is VERY early in development, but I'm sure they'd love  
 contributions!
 

Hi David,
That does look quite close to what I imagined, probably a good start then!
Romain


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] maximum value and corresponding index

2009-08-06 Thread Robert Kern
2009/8/6 Hans Meine me...@informatik.uni-hamburg.de:
 On Wednesday 05 August 2009 22:06:03 David Goldsmith wrote:
 But you can cheat and put them on one line (if that's all you're after):
  x = np.array([1, 2, 3])
  maxi = x.argmax(); maxv = x[maxi]

 Is there any reason not to put this as a convenience function into numpy?
 It is needed so frequently, and it's a shame that the shortest solution
 traverses the array twice (the above is usually longer due to variable names,
 and/or 1 dimensions).

The array is only traversed once. Indexing is O(1).

I'm -1 on adding a function to do this. If you want to keep that
convenience function in your own utilities, great.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Keith Goodman
On Thu, Aug 6, 2009 at 8:55 AM, josef.p...@gmail.com wrote:
 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them anymore?

 Josef

 a
 array([[1, 2, 3, 3, 0],
       [2, 2, 3, 2, 1]])
 a-a.max(0)
 array([[-1,  0,  0,  0, -1],
       [ 0,  0,  0, -1,  0]])
 a-a.max(1)
 Traceback (most recent call last):
  File pyshell#135, line 1, in module
    a-a.max(1)
 ValueError: shape mismatch: objects cannot be broadcast to a single shape
 a-a.max(1)[:,None]
 array([[-2, -1,  0,  0, -3],
       [-1, -1,  0, -1, -2]])

Would this do it?

 pylab.demean??
Type:   function
Base Class: type 'function'
String Form:function demean at 0x3c5c050
Namespace:  Interactive
File:   /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
Definition: pylab.demean(x, axis=0)
Source:
def demean(x, axis=0):
Return x minus its mean along the specified axis
x = np.asarray(x)
if axis:
ind = [slice(None)] * axis
ind.append(np.newaxis)
return x - x.mean(axis)[ind]
return x - x.mean(axis)
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] maximum value and corresponding index

2009-08-06 Thread Hans Meine
On Thursday 06 August 2009 17:27:51 Robert Kern wrote:
 2009/8/6 Hans Meine me...@informatik.uni-hamburg.de:
  On Wednesday 05 August 2009 22:06:03 David Goldsmith wrote:
  But you can cheat and put them on one line (if that's all you're 
after):
   x = np.array([1, 2, 3])
   maxi = x.argmax(); maxv = x[maxi]
 
  Is there any reason not to put this as a convenience function into numpy?
  It is needed so frequently, and it's a shame that the shortest solution
  traverses the array twice [...]

 The array is only traversed once. Indexing is O(1).

Yes, but I wrote the shortest solution, which is to call .argmax() and 
.max(), since..

  the above is usually longer due to variable names, and/or 1 dimensions

which requires unravel_index() (note that flattening is also inefficient when 
the array is strided).

Never mind,
  Hans
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 8:55 AM, josef.p...@gmail.com wrote:
 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them anymore?

 Josef

 a
 array([[1, 2, 3, 3, 0],
       [2, 2, 3, 2, 1]])
 a-a.max(0)
 array([[-1,  0,  0,  0, -1],
       [ 0,  0,  0, -1,  0]])
 a-a.max(1)
 Traceback (most recent call last):
  File pyshell#135, line 1, in module
    a-a.max(1)
 ValueError: shape mismatch: objects cannot be broadcast to a single shape
 a-a.max(1)[:,None]
 array([[-2, -1,  0,  0, -3],
       [-1, -1,  0, -1, -2]])

 Would this do it?

 pylab.demean??
 Type:           function
 Base Class:     type 'function'
 String Form:    function demean at 0x3c5c050
 Namespace:      Interactive
 File:           /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
 Definition:     pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
    Return x minus its mean along the specified axis
    x = np.asarray(x)
    if axis:
        ind = [slice(None)] * axis
        ind.append(np.newaxis)
        return x - x.mean(axis)[ind]
    return x - x.mean(axis)

Ouch! That doesn't handle axis=-1.

if axis != 0:
ind = [slice(None)] * x.ndim
ind[axis] = np.newaxis

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 11:15, Keith Goodmankwgood...@gmail.com wrote:
 thanksOn Thu, Aug 6, 2009 at 9:07 AM, Robert Kernrobert.k...@gmail.com 
 wrote:
 On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com wrote:

 pylab.demean??
 Type:           function
 Base Class:     type 'function'
 String Form:    function demean at 0x3c5c050
 Namespace:      Interactive
 File:           /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
 Definition:     pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
    Return x minus its mean along the specified axis
    x = np.asarray(x)
    if axis:
        ind = [slice(None)] * axis
        ind.append(np.newaxis)
        return x - x.mean(axis)[ind]
    return x - x.mean(axis)

 Ouch! That doesn't handle axis=-1.

 if axis != 0:
    ind = [slice(None)] * x.ndim
    ind[axis] = np.newaxis

 Hey, didn't you warn us about the dangers of if arr the other day?

Yes, but actually that wasn't quite the problem. if axis: would have
been fine, if a bit obscure, as long as the body was fixed to not
expect axis0.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread josef . pktd
On Thu, Aug 6, 2009 at 12:07 PM, Robert Kernrobert.k...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 8:55 AM, josef.p...@gmail.com wrote:
 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them anymore?

 Josef

 a
 array([[1, 2, 3, 3, 0],
       [2, 2, 3, 2, 1]])
 a-a.max(0)
 array([[-1,  0,  0,  0, -1],
       [ 0,  0,  0, -1,  0]])
 a-a.max(1)
 Traceback (most recent call last):
  File pyshell#135, line 1, in module
    a-a.max(1)
 ValueError: shape mismatch: objects cannot be broadcast to a single shape
 a-a.max(1)[:,None]
 array([[-2, -1,  0,  0, -3],
       [-1, -1,  0, -1, -2]])

 Would this do it?

 pylab.demean??
 Type:           function
 Base Class:     type 'function'
 String Form:    function demean at 0x3c5c050
 Namespace:      Interactive
 File:           /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
 Definition:     pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
    Return x minus its mean along the specified axis
    x = np.asarray(x)
    if axis:
        ind = [slice(None)] * axis
        ind.append(np.newaxis)
        return x - x.mean(axis)[ind]
    return x - x.mean(axis)

 Ouch! That doesn't handle axis=-1.

 if axis != 0:
    ind = [slice(None)] * x.ndim
    ind[axis] = np.newaxis


Thanks, that's it.

I have seen implementation of helper functions similar to this in
other packages, but I thought there is already something in numpy.  I
think this should be a simple helper function in numpy to avoid
mistakes and complicated implementation like the one in stats.nanstd
even if it's only a few lines.

Josef


 --
 Robert Kern

 I have come to believe that the whole world is an enigma, a harmless
 enigma that is made terrible by our own mad attempt to interpret it as
 though it had an underlying truth.
  -- Umberto Eco
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Keith Goodman
On Thu, Aug 6, 2009 at 9:18 AM, Robert Kernrobert.k...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 11:15, Keith Goodmankwgood...@gmail.com wrote:
 thanksOn Thu, Aug 6, 2009 at 9:07 AM, Robert Kernrobert.k...@gmail.com 
 wrote:
 On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com wrote:

 pylab.demean??
 Type:           function
 Base Class:     type 'function'
 String Form:    function demean at 0x3c5c050
 Namespace:      Interactive
 File:           /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
 Definition:     pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
    Return x minus its mean along the specified axis
    x = np.asarray(x)
    if axis:
        ind = [slice(None)] * axis
        ind.append(np.newaxis)
        return x - x.mean(axis)[ind]
    return x - x.mean(axis)

 Ouch! That doesn't handle axis=-1.

 if axis != 0:
    ind = [slice(None)] * x.ndim
    ind[axis] = np.newaxis

 Hey, didn't you warn us about the dangers of if arr the other day?

 Yes, but actually that wasn't quite the problem. if axis: would have
 been fine, if a bit obscure, as long as the body was fixed to not
 expect axis0.

Oh, of course. Thanks for pointing that out. Now I can fix a bug in my own code.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Pierre GM

On Aug 6, 2009, at 12:22 PM, Keith Goodman wrote:

 On Thu, Aug 6, 2009 at 9:18 AM, Robert Kernrobert.k...@gmail.com  
 wrote:
 On Thu, Aug 6, 2009 at 11:15, Keith Goodmankwgood...@gmail.com  
 wrote:
 thanksOn Thu, Aug 6, 2009 at 9:07 AM, Robert Kernrobert.k...@gmail.com 
  wrote:
 On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com  
 wrote:

 pylab.demean??
 Type:   function
 Base Class: type 'function'
 String Form:function demean at 0x3c5c050
 Namespace:  Interactive
 File:   /usr/lib/python2.6/dist-packages/matplotlib/ 
 mlab.py
 Definition: pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
Return x minus its mean along the specified axis
x = np.asarray(x)
if axis:
ind = [slice(None)] * axis
ind.append(np.newaxis)
return x - x.mean(axis)[ind]
return x - x.mean(axis)

FYI, there's a anom method for MaskedArrays that does the same thing  
as demean (if you can't /don't want to import mpl)
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 11:21, josef.p...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 12:07 PM, Robert Kernrobert.k...@gmail.com wrote:

 if axis != 0:
    ind = [slice(None)] * x.ndim
    ind[axis] = np.newaxis


 Thanks, that's it.

 I have seen implementation of helper functions similar to this in
 other packages, but I thought there is already something in numpy.  I
 think this should be a simple helper function in numpy to avoid
 mistakes and complicated implementation like the one in stats.nanstd
 even if it's only a few lines.

It would make a good contribution to numpy.lib.index_tricks.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Keith Goodman
thanksOn Thu, Aug 6, 2009 at 9:07 AM, Robert Kernrobert.k...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 11:03, Keith Goodmankwgood...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 8:55 AM, josef.p...@gmail.com wrote:
 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them anymore?

 Josef

 a
 array([[1, 2, 3, 3, 0],
       [2, 2, 3, 2, 1]])
 a-a.max(0)
 array([[-1,  0,  0,  0, -1],
       [ 0,  0,  0, -1,  0]])
 a-a.max(1)
 Traceback (most recent call last):
  File pyshell#135, line 1, in module
    a-a.max(1)
 ValueError: shape mismatch: objects cannot be broadcast to a single shape
 a-a.max(1)[:,None]
 array([[-2, -1,  0,  0, -3],
       [-1, -1,  0, -1, -2]])

 Would this do it?

 pylab.demean??
 Type:           function
 Base Class:     type 'function'
 String Form:    function demean at 0x3c5c050
 Namespace:      Interactive
 File:           /usr/lib/python2.6/dist-packages/matplotlib/mlab.py
 Definition:     pylab.demean(x, axis=0)
 Source:
 def demean(x, axis=0):
    Return x minus its mean along the specified axis
    x = np.asarray(x)
    if axis:
        ind = [slice(None)] * axis
        ind.append(np.newaxis)
        return x - x.mean(axis)[ind]
    return x - x.mean(axis)

 Ouch! That doesn't handle axis=-1.

 if axis != 0:
    ind = [slice(None)] * x.ndim
    ind[axis] = np.newaxis

Hey, didn't you warn us about the dangers of if arr the other day?
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Charles R Harris
On Thu, Aug 6, 2009 at 9:55 AM, josef.p...@gmail.com wrote:

 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them anymore?


Adding a keyword to retain the number of dimensions has been mooted. It
shouldn't be too difficult to implement and would allow things like:

 scaled = a/a.max(1, reduce=0)

I could do that for 1.4 if folks are interested.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread Keith Goodman
On Thu, Aug 6, 2009 at 9:58 AM, Charles R
Harrischarlesr.har...@gmail.com wrote:


 On Thu, Aug 6, 2009 at 9:55 AM, josef.p...@gmail.com wrote:

 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them
 anymore?

 Adding a keyword to retain the number of dimensions has been mooted. It
 shouldn't be too difficult to implement and would allow things like:

 scaled = a/a.max(1, reduce=0)

 I could do that for 1.4 if folks are interested.

I'd use that. It's better than what I usually do:

scaled = a / a.max(1).reshape(-1,1)
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread James Bergstra
David Warde-Farley dwf at cs.toronto.edu writes:
 It did inspire some of our colleagues in Montreal to create this,
 though:

      http://code.google.com/p/cuda-ndarray/

 I gather it is VERY early in development, but I'm sure they'd love
 contributions!


Hi David,
That does look quite close to what I imagined, probably a good start then!
Romain

Hi, I'm one of the devs for that project.   Thanks David for the link.
 I put some text on the homepage so it's a little more
self-explanatory.  We do welcome contributions.

I feel like I must be reinventing the wheel on this, so I'd really
appreciate it if someone who knows of a similar project would let me
know about it. Otherwise we'll keep plugging away at replicating core
ndarray interface elements (operators, math.h-type functions, array
indexing, etc.)

http://code.google.com/p/cuda-ndarray/

James
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] add axis to results of reduction (mean, min, ...)

2009-08-06 Thread josef . pktd
On Thu, Aug 6, 2009 at 1:07 PM, Keith Goodmankwgood...@gmail.com wrote:
 On Thu, Aug 6, 2009 at 9:58 AM, Charles R
 Harrischarlesr.har...@gmail.com wrote:


 On Thu, Aug 6, 2009 at 9:55 AM, josef.p...@gmail.com wrote:

 What's the best way of getting back the correct shape to be able to
 broadcast, mean, min,.. to the original array, that works for
 arbitrary dimension and axis?

 I thought I have seen some helper functions, but I don't find them
 anymore?

 Adding a keyword to retain the number of dimensions has been mooted. It
 shouldn't be too difficult to implement and would allow things like:

 scaled = a/a.max(1, reduce=0)

 I could do that for 1.4 if folks are interested.

 I'd use that. It's better than what I usually do:

 scaled = a / a.max(1).reshape(-1,1)
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


added feature to numpy reduce function would be nice.
but helper function is still useful for our own reduce operations


something like this function with an awful name

Josef



import numpy as np

def addreducedaxis(x,axis=None):
'''adds axis so that results of reduce operation broadcast to original
array


Parameter
-
x : array n-dim
   array that is the result of reduce operation, e.g. mean, min
axis : int
   axis that was removed in the reduce operation

Return
--

y : array (n+1)dim
   array with additional axis

'''
if axis != 0 and not axis is None:
ind = [slice(None)] * (x.ndim+1)
ind[axis] = np.newaxis
return x[ind]
else:
return x


a = np.array([[1,2,3,3,0],[2,2,3,2,1]])
a3 = np.dstack((a,a))

print np.all((a3-a3.mean(1)[:,None,:]) == (a3 - addreducedaxis(a3.mean(1),1)))
print np.all((a3-a3.mean(-2)[:,None,:]) == (a3 -
addreducedaxis(a3.mean(-2),-2)))
print np.all((a3-a3.mean()) == (a3 - addreducedaxis(a3.mean(

print np.all((a3.ravel()-a3.mean(None)) == (a3.ravel() -
addreducedaxis(a3.mean(None),None)))
print np.all((a3-a3.mean(None)) == (a3 - addreducedaxis(a3.mean(None),None)))


#example usage

from numpy.testing import assert_almost_equal

for axis in [None,0,1,2]:
m = a3.mean(axis)
v = ((a3 - addreducedaxis(m,axis))**2).mean(axis)
assert_almost_equal(v, np.var(a3,axis),15)

#normalize array along one axis
a3n = (a3 - 
addreducedaxis(np.mean(a3,1),1))/np.sqrt(addreducedaxis(np.var(a3,1),1))
print a3n.mean(1)
print a3n.var(1)
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Charles R Harris
On Thu, Aug 6, 2009 at 11:12 AM, James Bergstra
bergs...@iro.umontreal.cawrote:

 David Warde-Farley dwf at cs.toronto.edu writes:
  It did inspire some of our colleagues in Montreal to create this,
  though:
 
   http://code.google.com/p/cuda-ndarray/
 
  I gather it is VERY early in development, but I'm sure they'd love
  contributions!
 
 
 Hi David,
 That does look quite close to what I imagined, probably a good start then!
 Romain

 Hi, I'm one of the devs for that project.   Thanks David for the link.
  I put some text on the homepage so it's a little more
 self-explanatory.  We do welcome contributions.

 I feel like I must be reinventing the wheel on this, so I'd really
 appreciate it if someone who knows of a similar project would let me
 know about it. Otherwise we'll keep plugging away at replicating core
 ndarray interface elements (operators, math.h-type functions, array
 indexing, etc.)

 http://code.google.com/p/cuda-ndarray/


I almost looks like you are reimplementing numpy, in c++ no less. Is there
any reason why you aren't working with a numpy branch and just adding
ufuncs? I'm also curious if you have thoughts about how to use the GPU
pipelines in parallel.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Bug or documentation omission with dtype parameter with numpy.mean

2009-08-06 Thread Bruce Southey
Hi,
Should numpy.mean() (and similar functions) maintain the original dtype 
or the dtype used for the dtype option?
I do understand the numerical issues involved but my question relates to 
whether this should be a bug or needs clarification in the documentation.

According to the help, the dtype parameter to numpy.mean() says:
 dtype : dtype, optional
 Type to use in computing the mean. For integer inputs, the default
 is float64; for floating point, inputs it is the same as the input
 dtype.

I interpret this to indicate the type used in the internal calculations 
and technically does not say what the output dtype should be. But the 
dtype option does change the output dtype as the simple example below 
shows.

With Python 2.6 and numpy '1.4.0.dev7282' on Linux 64-bit Fedora 11

  import numpy as np
  a=np.array([1,2,3])
  a.dtype
dtype('int64')
  a.mean().dtype
dtype('float64')
  a.mean(dtype=np.float32).dtype
dtype('float64')
  a.mean(dtype=np.float64).dtype
dtype('float64')
  a.mean(dtype=np.float128).dtype
dtype('float128')
  a.mean(dtype=np.int).dtype
dtype('float64')

Clearly the output dtype is float64 or higher as determined by the input 
dtype or dtype parameter.

Bruce
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread James Bergstra
On Thu, Aug 6, 2009 at 1:19 PM, Charles R
Harrischarlesr.har...@gmail.com wrote:
 I almost looks like you are reimplementing numpy, in c++ no less. Is there
 any reason why you aren't working with a numpy branch and just adding
 ufuncs?

I don't know how that would work.  The Ufuncs need a datatype to work
with, and AFAIK, it would break everything if a numpy ndarray pointed
to memory on the GPU.  Could you explain what you mean a little more?

 I'm also curious if you have thoughts about how to use the GPU
 pipelines in parallel.

Current thinking for ufunc type computations:
1) divide up the tensors into subtensors whose dimensions have
power-of-two sizes (this permits a fast integer - ndarray coordinate
computation using bit shifting),
2) launch a kernel for each subtensor in it's own stream to use
parallel pipelines.
3) sync and return.

This is a pain to do without automatic code generation though.
Currently we're using macros, but that's not pretty.
C++ has templates, which we don't really use yet, but were planning on
using.  These have some power to generate code.
The 'theano' project (www.pylearn.org/theano) for which cuda-ndarray
was created has a more powerful code generation mechanism similar to
weave.   This algorithm is used in theano-cuda-ndarray.
Scipy.weave could be very useful for generating code for specific
shapes/ndims on demand, if weave could use nvcc.

James
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Erik Tollerud
Note that this is from a user perspective, as I have no particular plan of
developing the details of this implementation, but I've thought for a long
time that GPU support could be great for numpy (I would also vote for OpenCL
support over cuda, although conceptually they seem quite similar)...
But  what exactly would the large-scale plan be?  One of the advantages of
GPGPUs is that they are particularly suited to rather complicated
paralellizable algorithms, and the numpy-level basic operations are just the
simple arithmatic operations.  So while I'd love to see it working, it's
unclear to me exactly how much is gained at the core numpy level, especially
given that it's limited to single-precision on most GPUs.

Now linear algebra or FFTs on a GPU would probably be a huge boon, I'll
admit - especially if it's in the form of a drop-in replacement for the
numpy or scipy versions.

By the way, I noticed no one mentioned the GPUArray class in pycuda (and it
looks like there's something similar in the pyopencl) - seems like that's
already done a fair amount of the work...
http://documen.tician.de/pycuda/array.html#pycuda.gpuarray.GPUArray



On Thu, Aug 6, 2009 at 10:41 AM, James Bergstra
bergs...@iro.umontreal.cawrote:

 On Thu, Aug 6, 2009 at 1:19 PM, Charles R
 Harrischarlesr.har...@gmail.com wrote:
  I almost looks like you are reimplementing numpy, in c++ no less. Is
 there
  any reason why you aren't working with a numpy branch and just adding
  ufuncs?

 I don't know how that would work.  The Ufuncs need a datatype to work
 with, and AFAIK, it would break everything if a numpy ndarray pointed
 to memory on the GPU.  Could you explain what you mean a little more?

  I'm also curious if you have thoughts about how to use the GPU
  pipelines in parallel.

 Current thinking for ufunc type computations:
 1) divide up the tensors into subtensors whose dimensions have
 power-of-two sizes (this permits a fast integer - ndarray coordinate
 computation using bit shifting),
 2) launch a kernel for each subtensor in it's own stream to use
 parallel pipelines.
 3) sync and return.

 This is a pain to do without automatic code generation though.
 Currently we're using macros, but that's not pretty.
 C++ has templates, which we don't really use yet, but were planning on
 using.  These have some power to generate code.
 The 'theano' project (www.pylearn.org/theano) for which cuda-ndarray
 was created has a more powerful code generation mechanism similar to
 weave.   This algorithm is used in theano-cuda-ndarray.
 Scipy.weave could be very useful for generating code for specific
 shapes/ndims on demand, if weave could use nvcc.

 James
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Matthieu Brucher
2009/8/6 Erik Tollerud erik.tolle...@gmail.com:
 Note that this is from a user perspective, as I have no particular plan of
 developing the details of this implementation, but I've thought for a long
 time that GPU support could be great for numpy (I would also vote for OpenCL
 support over cuda, although conceptually they seem quite similar)...
 But  what exactly would the large-scale plan be?  One of the advantages of
 GPGPUs is that they are particularly suited to rather complicated
 paralellizable algorithms,

You mean simple parallizable algorithms, I suppose?

 and the numpy-level basic operations are just the
 simple arithmatic operations.  So while I'd love to see it working, it's
 unclear to me exactly how much is gained at the core numpy level, especially
 given that it's limited to single-precision on most GPUs.
 Now linear algebra or FFTs on a GPU would probably be a huge boon, I'll
 admit - especially if it's in the form of a drop-in replacement for the
 numpy or scipy versions.
 By the way, I noticed no one mentioned the GPUArray class in pycuda (and it
 looks like there's something similar in the pyopencl) - seems like that's
 already done a fair amount of the work...
 http://documen.tician.de/pycuda/array.html#pycuda.gpuarray.GPUArray


 On Thu, Aug 6, 2009 at 10:41 AM, James Bergstra bergs...@iro.umontreal.ca
 wrote:

 On Thu, Aug 6, 2009 at 1:19 PM, Charles R
 Harrischarlesr.har...@gmail.com wrote:
  I almost looks like you are reimplementing numpy, in c++ no less. Is
  there
  any reason why you aren't working with a numpy branch and just adding
  ufuncs?

 I don't know how that would work.  The Ufuncs need a datatype to work
 with, and AFAIK, it would break everything if a numpy ndarray pointed
 to memory on the GPU.  Could you explain what you mean a little more?

  I'm also curious if you have thoughts about how to use the GPU
  pipelines in parallel.

 Current thinking for ufunc type computations:
 1) divide up the tensors into subtensors whose dimensions have
 power-of-two sizes (this permits a fast integer - ndarray coordinate
 computation using bit shifting),
 2) launch a kernel for each subtensor in it's own stream to use
 parallel pipelines.
 3) sync and return.

 This is a pain to do without automatic code generation though.
 Currently we're using macros, but that's not pretty.
 C++ has templates, which we don't really use yet, but were planning on
 using.  These have some power to generate code.
 The 'theano' project (www.pylearn.org/theano) for which cuda-ndarray
 was created has a more powerful code generation mechanism similar to
 weave.   This algorithm is used in theano-cuda-ndarray.
 Scipy.weave could be very useful for generating code for specific
 shapes/ndims on demand, if weave could use nvcc.

 James
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion



 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion





-- 
Information System Engineer, Ph.D.
Website: http://matthieu-brucher.developpez.com/
Blogs: http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn: http://www.linkedin.com/in/matthieubrucher
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Optimized half-sizing of images?

2009-08-06 Thread Christopher Barker
(second try on this message. the fist time I included a test PNG that 
made it too large)


Hi folks,

We have a need to to generate half-size version of RGB images as quickly
as possible. PIL does a pretty good job, but it dawned on me that in the
special case of a half-size, one might be able to do it faster with
numpy, simply averaging the four pixels in the larger image to create
one in the small. I'm doing tiling, and thus reducing 512x512 images to
256x256, so I imagine I'm making good use of cache (it does get pretty
pokey with really large images!)


What I have now is essentially this :

# a is a (h, w, 3) RGB array

 a2 =  a[0::2, 0::2, :].astype(np.uint16)
 a2 += a[0::2, 1::2, :]
 a2 += a[1::2, 0::2, :]
 a2 += a[1::2, 1::2, :]
 a2 /= 4
 return a2.astype(np.uint8)

time:  67.2 ms per loop

I can speed it up a bit if I accumulate in a uint8 and divide as I go to
prevent overflow:

 a2 =  a[0::2, 0::2, :].astype(np.uint8) / 4
 a2 += a[0::2, 1::2, :] / 4
 a2 += a[1::2, 0::2, :] / 4
 a2 += a[1::2, 1::2, :] / 4
 return a2

time:  46.6 ms per loop

That does lose a touch of accuracy, I suppose, but nothing I can see.


Method 1 is about twice as slow as PIL's bilinear scaling.

Can I do better? It seems it should be faster if I can avoid so many
separate loops through the array. I figure there may be some way with
filter or convolve or ndimage, but they all seem to return an array the
same size.

Any ideas?

(Cython is another option, of course)

-Chris

Test code  enclosed.









--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/ORR(206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115   (206) 526-6317   main reception

chris.bar...@noaa.gov



numpy_resize.py
Description: application/python
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Optimized half-sizing of images?

2009-08-06 Thread Stéfan van der Walt
Hi Chris

2009/8/6 Christopher Barker chris.bar...@noaa.gov:
 Can I do better? It seems it should be faster if I can avoid so many
 separate loops through the array. I figure there may be some way with
 filter or convolve or ndimage, but they all seem to return an array the
 same size.

Are you willing to depend on SciPy?  We've got pretty fast zooming
code in ndimage.

If speed is a big issue, I'd consider using the GPU, which was made
for this sort of down-sampling.

Cheers
Stéfan
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread David Warde-Farley
On 6-Aug-09, at 2:54 PM, Erik Tollerud wrote:

 Now linear algebra or FFTs on a GPU would probably be a huge boon,  
 I'll
 admit - especially if it's in the form of a drop-in replacement for  
 the
 numpy or scipy versions.

The word I'm hearing from people in my direct acquaintance who are  
using it is that if you have code that even do lots of matrix  
multiplies, nevermind solving systems or anything like that, the  
speedup is several orders of magnitude. Things that used to take weeks  
now take a day or two. If you can deal with the loss of precision it's  
really quite worth it.

 By the way, I noticed no one mentioned the GPUArray class in pycuda  
 (and it
 looks like there's something similar in the pyopencl) - seems like  
 that's
 already done a fair amount of the work...
 http://documen.tician.de/pycuda/array.html#pycuda.gpuarray.GPUArray


This seems like a great start, I agree. The lack of any documentation  
on 'dot' is worrying, though.

David
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden

 Now linear algebra or FFTs on a GPU would probably be a huge boon, 
 I'll admit - especially if it's in the form of a drop-in replacement 
 for the numpy or scipy versions.


NumPy generate temporary arrays for expressions involving ndarrays. This 
extra allocation and copying often takes more time than the computation. 
With GPGPUs, we have to bus the data to and from VRAM as well. D. Knuth 
quoted Hoare saying that premature optimization is the root of all 
evil. Optimizing computation when the bottleneck is memory is premature.

In order to improve on this, I think we have to add lazy evaluation to 
NumPy. That is, an operator should not return a temporary array but a 
symbolic expression. So if we have an expression like

y = a*x + b

it should not evalute a*x into a temporary array. Rather, the operators 
would build up a parse tree like

y = add(multiply(a,x),b)

and evalute the whole expression  later on.

This would require two things: First we need dynamic code generation, 
which incidentally is what OpenCL is all about. I.e. OpenCL is 
dynamically invoked compiler; there is a function 
clCreateProgramFromSource, which  does just what it says. Second, we 
need arrays to be immutable. This is very important. If arrays are not 
immutable, code like this could fail:

y = a*x + b
x[0] = 1235512371235

With lazy evaluation, the memory overhead would be much smaller. The 
GPGPU would also get a more complex expressions to use as a kernels.

There should be an option of running this on the CPU, possibly using 
OpenMP for multi-threading. We could either depend on a compiler (C or 
Fortran) being installed, or use opcodes for a dedicated virtual machine 
(cf. what numexpr does).

In order to reduce the effect of immutable arrays, we could introduce a 
context-manager. Inside the with statement, all arrays would be 
immutable. Second, the __exit__ method could trigger the code generator 
and do all the evaluation. So we would get something like this:

# normal numpy here

with numpy.accelerator():

# arrays become immutable
# lazy evaluation
   
# code generation and evaluation on exit

# normal numpy continues here


Thus, here is my plan:

1. a special context-manager class
2. immutable arrays inside with statement
3. lazy evaluation: expressions build up a parse tree
4. dynamic code generation
5. evaluation on exit

I guess it is possibly to find ways to speed up this as well. If a 
context manager would always generate the same OpenCL code, the with 
statement would only need to execute once (we could raise an exception 
on enter to jump directly to exit).

It is possibly to create a superfast NumPy. But just plugging GPGPUs 
into the current design would be premature. In NumPy's current state, 
with mutable ndarrays and operators generating temporary arrays, there 
is not much to gain from introducing GPGPUs. It would only be beneficial 
in computationally demanding parts like FFTs and solvers for linear 
algebra and differential equations. Ufuncs with trancendental functions 
might also benefit. SciPy would certainly benefit more from GPGPUs than 
NumPy.

Just my five cents :-)

Regards,
Sturla Molden
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Optimized half-sizing of images?

2009-08-06 Thread Christopher Barker
Stéfan van der Walt wrote:
 Are you willing to depend on SciPy?  We've got pretty fast zooming
 code in ndimage.

I looked there, and didn't see it -- didn't think to look for zoom. 
Now I do, I'll give it a try.

However, my thought is that for the special case of half-sizing, all 
that spline stuff could be unneeded and slower. I'll see what we get.

thanks,

-Chris

-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/ORR(206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115   (206) 526-6317   main reception

chris.bar...@noaa.gov
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 15:57, Sturla Moldenstu...@molden.no wrote:

 Now linear algebra or FFTs on a GPU would probably be a huge boon,
 I'll admit - especially if it's in the form of a drop-in replacement
 for the numpy or scipy versions.

 NumPy generate temporary arrays for expressions involving ndarrays. This
 extra allocation and copying often takes more time than the computation.
 With GPGPUs, we have to bus the data to and from VRAM as well. D. Knuth
 quoted Hoare saying that premature optimization is the root of all
 evil. Optimizing computation when the bottleneck is memory is premature.

 It is possibly to create a superfast NumPy. But just plugging GPGPUs
 into the current design would be premature. In NumPy's current state,
 with mutable ndarrays and operators generating temporary arrays, there
 is not much to gain from introducing GPGPUs. It would only be beneficial
 in computationally demanding parts like FFTs and solvers for linear
 algebra and differential equations.

I believe that is exactly the point that Erik is making. :-)

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden
Robert Kern wrote:
 I believe that is exactly the point that Erik is making. :-)
   
I wasn't arguing against him, just suggesting a solution. :-)

I have big hopes for lazy evaluation, if we can find a way to to it right.

Sturla
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread James Bergstra
On Thu, Aug 6, 2009 at 4:57 PM, Sturla Moldenstu...@molden.no wrote:

 Now linear algebra or FFTs on a GPU would probably be a huge boon,
 I'll admit - especially if it's in the form of a drop-in replacement
 for the numpy or scipy versions.


 NumPy generate temporary arrays for expressions involving ndarrays. This
 extra allocation and copying often takes more time than the computation.
 With GPGPUs, we have to bus the data to and from VRAM as well. D. Knuth
 quoted Hoare saying that premature optimization is the root of all
 evil. Optimizing computation when the bottleneck is memory is premature.

 In order to improve on this, I think we have to add lazy evaluation to
 NumPy. That is, an operator should not return a temporary array but a
 symbolic expression. So if we have an expression like

    y = a*x + b

 it should not evalute a*x into a temporary array. Rather, the operators
 would build up a parse tree like

    y = add(multiply(a,x),b)

 and evalute the whole expression  later on.
[snip]
 Regards,
 Sturla Molden
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


Hi Sturla,

The plan you describe is a good one, and Theano
(www.pylearn.org/theano) almost exactly implements it.  You should
check it out.  It does not use 'with' syntax at the moment, but it
could provide the backend machinery for your mechanism if you want to
go forward with that.  Theano provides
- symbolic expression building for a big subset of what numpy can do
(and a few things that it doesn't)
- expression optimization (for faster and more accurate computations)
- dynamic code generation
- cacheing of compiled functions to disk.

Also, when you have a symbolic expression graph you can do cute stuff
like automatic differentiation.  We're currently working on the bridge
between theano and cuda so that you declare certain inputs as residing
on the GPU instead of the host memory, so you don't have to transfer
things to and from host memory as much.

James
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] FCompiler and runtime library dirs

2009-08-06 Thread Lisandro Dalcin
Hi, folks, Using NumPy 1.3.0 from Fedora 11, though this issue likely
applies to current trunk (I've not actually tested, just taken a look
at the sources)

As numpy.distutils.FCompiler inherits from
distutils.ccompiler.CCompiler, the method
runtime_library_dir_option() fails with NotImplementedError. I had
to add the monkeypatch pasted below to a setup.py script (full code at
 http://petsc.cs.iit.edu/petsc4py/petsc4py-dev/file/tip/demo/wrap-f2py/setup.py)
in order to get things working (wiht GCC on Linux):

from numpy.distutils.fcompiler import FCompiler
from numpy.distutils.unixccompiler import UnixCCompiler
FCompiler.runtime_library_dir_option = \
UnixCCompiler.runtime_library_dir_option.im_func

Do any of you have an idea about how to properly fix this issue in
numpy? I'm tempted to re-use the UnixCompiler implementation in POSIX
(including Mac OS X?), just to save some lines and do not repeat that
code in every FCompiler subclass... Comments?


-- 
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
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Strange Error with NumPy

2009-08-06 Thread Nanime Puloski
Can anyone explain to me why I receive an error(AtrributeError) in NumPy
when I do numpy.sin(2**64), but not when I do numpy.sin(2.0**64),
numpy.sin(float(2**64)) or even numpy.sin(2)?

Where is the root problem in all of this?
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Charles R Harris
On Thu, Aug 6, 2009 at 3:29 PM, James Bergstra bergs...@iro.umontreal.cawrote:

 On Thu, Aug 6, 2009 at 4:57 PM, Sturla Moldenstu...@molden.no wrote:
 
  Now linear algebra or FFTs on a GPU would probably be a huge boon,
  I'll admit - especially if it's in the form of a drop-in replacement
  for the numpy or scipy versions.
 
 
  NumPy generate temporary arrays for expressions involving ndarrays. This
  extra allocation and copying often takes more time than the computation.
  With GPGPUs, we have to bus the data to and from VRAM as well. D. Knuth
  quoted Hoare saying that premature optimization is the root of all
  evil. Optimizing computation when the bottleneck is memory is premature.
 
  In order to improve on this, I think we have to add lazy evaluation to
  NumPy. That is, an operator should not return a temporary array but a
  symbolic expression. So if we have an expression like
 
 y = a*x + b
 
  it should not evalute a*x into a temporary array. Rather, the operators
  would build up a parse tree like
 
 y = add(multiply(a,x),b)
 
  and evalute the whole expression  later on.
 [snip]
  Regards,
  Sturla Molden
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion
 

 Hi Sturla,

 The plan you describe is a good one, and Theano
 (www.pylearn.org/theano) almost exactly implements it.  You should
 check it out.  It does not use 'with' syntax at the moment, but it
 could provide the backend machinery for your mechanism if you want to
 go forward with that.  Theano provides
 - symbolic expression building for a big subset of what numpy can do
 (and a few things that it doesn't)
 - expression optimization (for faster and more accurate computations)
 - dynamic code generation
 - cacheing of compiled functions to disk.

 Also, when you have a symbolic expression graph you can do cute stuff
 like automatic differentiation.  We're currently working on the bridge
 between theano and cuda so that you declare certain inputs as residing
 on the GPU instead of the host memory, so you don't have to transfer
 things to and from host memory as much.


So what simple things could numpy implement that would help here? It almost
sounds like numpy would mostly be an interface to python and the gpu would
execute specialized code written and compiled for specific problems. Whether
the code that gets compiled is written using lazy evaluation (ala Sturla),
or is expressed some other way seems like an independent issue. It sounds
like one important thing would be having arrays that reside on the GPU.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Strange Error with NumPy

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 17:01, Nanime Puloskinpulo...@gmail.com wrote:
 Can anyone explain to me why I receive an error(AtrributeError) in NumPy
 when I do numpy.sin(2**64), but not when I do numpy.sin(2.0**64),
 numpy.sin(float(2**64)) or even numpy.sin(2)?

 Where is the root problem in all of this?

numpy deals with objects that can be natively represented by the C
numerical types on your machine. On your machine 2**64 gives you a
Python long object which cannot be converted to one of the native C
numerical types, so numpy.sin() treats it as a regular Python object.
The default implementation of all of the ufuncs when faced with a
Python object it doesn't know about is to look for a method on the
object of the same name. long.sin() does not exist, so you get an
AttributeError.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Strange Error with NumPy

2009-08-06 Thread David Warde-Farley
On 6-Aug-09, at 6:01 PM, Nanime Puloski wrote:

 Can anyone explain to me why I receive an error(AtrributeError) in  
 NumPy
 when I do numpy.sin(2**64), but not when I do numpy.sin(2.0**64),
 numpy.sin(float(2**64)) or even numpy.sin(2)?

In [6]: type(2**64)
Out[6]: type 'long'

In [7]: type(2)
Out[7]: type 'int'

In [8]: type(2.0**64)
Out[8]: type 'float'

Probably because 2**64 yields a long, which is an arbitrary precision  
type in Python, and numpy is having trouble casting it to something it  
can use (it's too big for numpy.int64). Notice that np.sin(2**63 - 1)  
works fine while np.sin(2**63) doesn't - 2**63 - 1 is the largest  
value that an int64 can hold.

You're right that it could be approximately represented as a float32  
or float64, but numpy won't cast from exact types to inexact types  
without you telling it to do so.

I agree that the error message is less than ideal.

David
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Optimized half-sizing of images?

2009-08-06 Thread Christopher Barker

Christopher Barker wrote:

Stéfan van der Walt wrote:

Are you willing to depend on SciPy?  We've got pretty fast zooming
code in ndimage.


However, my thought is that for the special case of half-sizing, all 
that spline stuff could be unneeded and slower. I'll see what we get.


I've given that a try. The docs are pretty sparse. I couldn't figure out 
a way to get it to do the whole image at once. Rather I had to do each 
band separately. It ended up about the same speed as my int accumulator 
method:


def test5(a):

using ndimage

tested on 512x512 RGB image:
  time:  40 ms per loop


h = a.shape[0]/2
w = a.shape[1]/2
a2 = np.empty((h, w, 3), dtype=np.uint8)
a2[:,:,0] = ndi.zoom(a[:,:,0], 0.5, order=1)
a2[:,:,1] = ndi.zoom(a[:,:,1], 0.5, order=1)
a2[:,:,2] = ndi.zoom(a[:,:,2], 0.5, order=1)
return a2

Is there a way to get it to do all three colorbands at once?

NOTE: if I didn't set order to 1, it was MUCH slower, as one might expect.






--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/ORR(206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115   (206) 526-6317   main reception

chris.bar...@noaa.gov


numpy_resize.py
Description: application/python
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden
Charles R Harris wrote:
 Whether the code that gets compiled is written using lazy evaluation 
 (ala Sturla), or is expressed some other way seems like an independent 
 issue. It sounds like one important thing would be having arrays that 
 reside on the GPU.
Memory management is slow compared to computation. Operations like 
malloc, free and memcpy is not faster for VRAM than for RAM. There will 
be no benefit from the GPU if the bottleneck is memory. That is why we 
need to get rid of the creation of temporary arrays, hence lazy evaluation.

Having arrays reside in VRAM would reduce the communication between RAM 
and VRAM, but the problem with temporary arrays is still there.

Also VRAM tends to be a limited resource.

Sturla

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden
Sturla Molden wrote:
 Memory management is slow compared to computation. Operations like 
 malloc, free and memcpy is not faster for VRAM than for RAM. 

Actually it's not VRAM anymore, but whatever you call the memory 
dedicated to the GPU.

It is cheap to put 8 GB of RAM into a computer, but graphics cards with 
more than 1 GB memory are expensive and uncommon on e.g. laptops. And 
this memory will be needed for other things as well, e.g. graphics.

Sturla



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Charles R Harris
On Thu, Aug 6, 2009 at 4:36 PM, Sturla Molden stu...@molden.no wrote:

 Charles R Harris wrote:
  Whether the code that gets compiled is written using lazy evaluation
  (ala Sturla), or is expressed some other way seems like an independent
  issue. It sounds like one important thing would be having arrays that
  reside on the GPU.
 Memory management is slow compared to computation. Operations like
 malloc, free and memcpy is not faster for VRAM than for RAM. There will
 be no benefit from the GPU if the bottleneck is memory. That is why we
 need to get rid of the creation of temporary arrays, hence lazy evaluation.

 Having arrays reside in VRAM would reduce the communication between RAM
 and VRAM, but the problem with temporary arrays is still there.


I'm not arguing with that, but I regard it as a separate problem. One could,
after all, simply use an expression to GPU compiler to generate modules. The
question is what simple additions we can make to numpy so that it acts as a
convenient io channel. I mean, once the computations are moved elsewhere
numpy is basically a convenient way to address memory.



 Also VRAM tends to be a limited resource.


But getting less so. These days it comes in gigabytes and there is no reason
why it shouldn't soon excede what many folks have for main memory.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden
Charles R Harris wrote:

 I mean, once the computations are moved elsewhere numpy is basically a 
 convenient way to address memory.

That is how I mostly use NumPy, though. Computations I often do in 
Fortran 95 or C.

NumPy arrays on the GPU memory is an easy task. But then I would have to 
write the computation in OpenCL's dialect of C99? But I'd rather program 
everything in Python if I could. Details like GPU and OpenCL should be 
hidden away. Nice looking Python with NumPy is much easier to read and 
write. That is why I'd like to see a code generator (i.e. JIT compiler) 
for NumPy.


Sturla


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] Strange Error with NumPy Addendum

2009-08-06 Thread Nanime Puloski
Thank you for your responses so far.
What I also do not understand is why sin(2**64) works with
the standard Python math module, but fails to do so with NumPy?

To Robert Kern:
Can't 2^64 be represented in C as a long double?
It seems to work well on my machine.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Sturla Molden
James Bergstra wrote:
 The plan you describe is a good one, and Theano
 (www.pylearn.org/theano) almost exactly implements it.  You should
 check it out.  It does not use 'with' syntax at the moment, but it
 could provide the backend machinery for your mechanism if you want to
 go forward with that.  Theano provides
 - symbolic expression building for a big subset of what numpy can do
 (and a few things that it doesn't)
 - expression optimization (for faster and more accurate computations)
 - dynamic code generation
 - cacheing of compiled functions to disk.
Thank you James, theano looks great. :-D

Sturla


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Charles R Harris
On Thu, Aug 6, 2009 at 5:10 PM, Sturla Molden stu...@molden.no wrote:

 Charles R Harris wrote:

  I mean, once the computations are moved elsewhere numpy is basically a
  convenient way to address memory.

 That is how I mostly use NumPy, though. Computations I often do in
 Fortran 95 or C.

 NumPy arrays on the GPU memory is an easy task.


Glad to hear it. So maybe some way to specify and track where the memory is
allocated would be helpful. Travis wants to add a dictionary to ndarrays and
that might be useful here.

But then I would have to
 write the computation in OpenCL's dialect of C99? But I'd rather program
 everything in Python if I could. Details like GPU and OpenCL should be
 hidden away. Nice looking Python with NumPy is much easier to read and
 write. That is why I'd like to see a code generator (i.e. JIT compiler)
 for NumPy.


Yes, but that is a language/compiler problem. I'm thinking of what tools
numpy can offer that would help people experimenting with different
approaches to using GPUs. At some point we might want to adopt a working
approach but now seems a bit early for that.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Strange Error with NumPy Addendum

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 18:26, Nanime Puloskinpulo...@gmail.com wrote:
 Thank you for your responses so far.
 What I also do not understand is why sin(2**64) works with
 the standard Python math module, but fails to do so with NumPy?

math.sin() always converts the argument to a float. We do not.

 To Robert Kern:
 Can't 2^64 be represented in C as a long double?

For that value, yes, but not for long objects in general. We don't
look at the value itself, just the type.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Fernando Perez
On Thu, Aug 6, 2009 at 1:57 PM, Sturla Moldenstu...@molden.no wrote:
 In order to reduce the effect of immutable arrays, we could introduce a
 context-manager. Inside the with statement, all arrays would be
 immutable. Second, the __exit__ method could trigger the code generator
 and do all the evaluation. So we would get something like this:

    # normal numpy here

    with numpy.accelerator():

        # arrays become immutable
        # lazy evaluation

        # code generation and evaluation on exit

    # normal numpy continues here


 Thus, here is my plan:

 1. a special context-manager class
 2. immutable arrays inside with statement
 3. lazy evaluation: expressions build up a parse tree
 4. dynamic code generation
 5. evaluation on exit

You will face one issue here: unless you raise a special exception
inside the with block, the python interpreter will unconditionally
execute that code without your control.  I had a long talk about this
with Alex Martelli last year at scipy, where I pitched the idea of
allowing context managers to have an optional third method,
__execute__, which would get the code block in the with statement for
execution.  He was fairly pessimistic about the possibility of this
making its way into python, mostly (if I recall correctly) because of
scoping issues: the with statement does not introduce a new scope, so
you'd need to pass to this method the code plus the locals/globals of
the entire enclosing scope, which felt messy.

There was also the thorny question of how to pass the code block.
Source? Bytecode? What?  In many environments the source may not be
available.  Last year I wrote a gross hack to do this, which you can
find here:

http://bazaar.launchpad.net/~ipython-dev/ipython/0.10/annotate/head%3A/IPython/kernel/contexts.py

The idea is that it would be used by code like this (note, this
doesn't actually work right now):

def test_simple():

# XXX - for now, we need a running cluster to be started separately.  The
# daemon work is almost finished, and will make much of this unnecessary.
from IPython.kernel import client
mec = client.MultiEngineClient(('127.0.0.1',10105))

try:
mec.get_ids()
except ConnectionRefusedError:
import os, time
os.system('ipcluster -n 2 ')
time.sleep(2)
mec = client.MultiEngineClient(('127.0.0.1',10105))

mec.block = False

parallel = RemoteMultiEngine(mec)

mec.pushAll()

with parallel as pr:
# A comment
remote()  # this means the code below only runs remotely
print 'Hello remote world'
x = range(10)
# Comments are OK
# Even misindented.
y = x+1

print pr.x + pr.y

###

The problem with my approach is that I find it brittle and ugly enough
that I ultimately abandoned it.  I'd love to see if you find a proper
solution for this...

Cheers,

f
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Optimized half-sizing of images?

2009-08-06 Thread Zachary Pincus
 We have a need to to generate half-size version of RGB images as  
 quickly
 as possible.

How good do these need to look? You could just throw away every other  
pixel... image[::2, ::2].

Failing that, you could also try using ndimage's convolve routines to  
run a 2x2 box filter over the image, and then throw away half of the  
pixels. But this would be slower than optimal, because the kernel  
would be convolved over every pixel, not just the ones you intend to  
keep.

Really though, I'd just bite the bullet and write a C extension (or  
cython, whatever, an extension to work for a defined-dimensionality,  
defined-dtype array is pretty simple), or as suggested before, do it  
on the GPU. (Though I find that readback from the GPU can be slow  
enough that C code can beat it in some cases.)

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Fwd: GPU Numpy

2009-08-06 Thread Robert Kern
On Thu, Aug 6, 2009 at 19:00, Fernando Perezfperez@gmail.com wrote:
 On Thu, Aug 6, 2009 at 1:57 PM, Sturla Moldenstu...@molden.no wrote:
 In order to reduce the effect of immutable arrays, we could introduce a
 context-manager. Inside the with statement, all arrays would be
 immutable. Second, the __exit__ method could trigger the code generator
 and do all the evaluation. So we would get something like this:

    # normal numpy here

    with numpy.accelerator():

        # arrays become immutable
        # lazy evaluation

        # code generation and evaluation on exit

    # normal numpy continues here


 Thus, here is my plan:

 1. a special context-manager class
 2. immutable arrays inside with statement
 3. lazy evaluation: expressions build up a parse tree
 4. dynamic code generation
 5. evaluation on exit

 You will face one issue here: unless you raise a special exception
 inside the with block, the python interpreter will unconditionally
 execute that code without your control.  I had a long talk about this
 with Alex Martelli last year at scipy, where I pitched the idea of
 allowing context managers to have an optional third method,
 __execute__, which would get the code block in the with statement for
 execution.  He was fairly pessimistic about the possibility of this
 making its way into python, mostly (if I recall correctly) because of
 scoping issues: the with statement does not introduce a new scope, so
 you'd need to pass to this method the code plus the locals/globals of
 the entire enclosing scope, which felt messy.

Sometimes, I fantasize about writing a python4ply grammar that
repurposes the `` quotes to provide expression literals and ``` ```
triple quotes for multiline statement literals. They would be literals
for _ast abstract syntax trees.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] datetime code

2009-08-06 Thread Charles R Harris
Travis, Robert,

Is there any reason not to merge the c code in the datetime branch at this
time? If not, I will do it.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Strange Error with NumPy Addendum

2009-08-06 Thread David Warde-Farley
On 6-Aug-09, at 7:29 PM, Robert Kern wrote:

 For that value, yes, but not for long objects in general. We don't
 look at the value itself, just the type.

Err, don't look at the value (of a long), except when it's  
representable with an integer dtype, right? Hence why 2**63 - 1 works.

David
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Strange Error with NumPy Addendum

2009-08-06 Thread Robert Kern
On Fri, Aug 7, 2009 at 00:03, David Warde-Farleyd...@cs.toronto.edu wrote:
 On 6-Aug-09, at 7:29 PM, Robert Kern wrote:

 For that value, yes, but not for long objects in general. We don't
 look at the value itself, just the type.

 Err, don't look at the value (of a long), except when it's
 representable with an integer dtype, right? Hence why 2**63 - 1 works.

Err, you may be right.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion