[Numpy-discussion] A correction to numpy trapz function

[Nadav Horesh]

Here is what I get with the orriginal trapz function:

IDLE 1.2.2      
>>> import numpy as np
>>> np.__version__
'1.1.0'
>>> y = np.arange(24).reshape(6,4)
>>> x = np.arange(6)
>>> np.trapz(y, x, axis=0)

Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    np.trapz(y, x, axis=0)
  File "C:\Python25\Lib\site-packages\numpy\lib\function_base.py", line 1536, in trapz
    return add.reduce(d * (y[slice1]+y[slice2])/2.0,axis)
ValueError: shape mismatch: objects cannot be broadcast to a single shape
>>> 

   Nadav.

-----הודעה מקורית-----
מאת: numpy-discussion-bounces@scipy.org בשם Ryan May
נשלח: ש 12-יולי-08 18:31
אל: Discussion of Numerical Python
נושא: Re: [Numpy-discussion] A correction to numpy trapz function
 
Nadav Horesh wrote:
> The function trapz accepts x axis vector only for axis=-1. Here is my modification (correction?) to let it accept a vector x for integration along any axis:
> 
> def trapz(y, x=None, dx=1.0, axis=-1):
>     """
>     Integrate y(x) using samples along the given axis and the composite
>     trapezoidal rule.  If x is None, spacing given by dx is assumed. If x
>     is an array, it must have either the dimensions of y, or a vector of
>     length matching the dimension of y along the integration axis.
>     """
>     y = asarray(y)
>     nd = y.ndim
>     slice1 = [slice(None)]*nd
>     slice2 = [slice(None)]*nd
>     slice1[axis] = slice(1,None)
>     slice2[axis] = slice(None,-1)
>     if x is None:
>         d = dx
>     else:
>         x = asarray(x)
>         if x.ndim == 1:
>             if len(x) != y.shape[axis]:
>                 raise ValueError('x length (%d) does not match y axis %d length (%d)' % (len(x), axis, y.shape[axis]))
>             d = diff(x)
>             return tensordot(d, (y[slice1]+y[slice2])/2.0,(0, axis))
>         d = diff(x, axis=axis)
>     return add.reduce(d * (y[slice1]+y[slice2])/2.0,axis)
> 

What version were you working with originally? With 1.1, this is what I 
have:

def trapz(y, x=None, dx=1.0, axis=-1):
     """Integrate y(x) using samples along the given axis and the composite
     trapezoidal rule.  If x is None, spacing given by dx is assumed.
     """
     y = asarray(y)
     if x is None:
         d = dx
     else:
         d = diff(x,axis=axis)
     nd = len(y.shape)
     slice1 = [slice(None)]*nd
     slice2 = [slice(None)]*nd
     slice1[axis] = slice(1,None)
     slice2[axis] = slice(None,-1)
     return add.reduce(d * (y[slice1]+y[slice2])/2.0,axis)

For me, this works fine with supplying x for axis != -1.

Ryan

-- 
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 4084 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/numpy-discussion/attachments/20080712/4dc16862/attachment.bin