A correction to numpy trapz function
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) Nadav.
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
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
Nadav Horesh wrote:
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
(Try not to top post on this list.) I can get it to work like this: import numpy as np y = np.arange(24).reshape(6,4) x = np.arange(6).reshape(-1,1) np.trapz(y, x, axis=0) From the text of the error message, you can see this is a problem with broadcasting. Due to broadcasting rules (which will *prepend* dimensions with size 1), you need to manually add an extra dimension to the end. Once I resize x, I can get this to work. You might want to look at this: http://www.scipy.org/EricsBroadcastingDoc Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma
I am aware that the error is related to the broadcasting, and that it can be solved by matching the shape of x to that of y --- this is how I solved it in the first place. I was thinking that the function "promises" to integrate over an array given a x vector and the axis, so let obscure the broadcasting rules and just enable it to do the work. There is a reason to leave trapz as it is (or even drop it) since numpy should stay as close as possible to "bare metal", but this function is borrowed also by scipy integration package, thus I rather give it a face lift. Nadav. -----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם Ryan May נשלח: ש 12-יולי-08 22:24 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] A correction to numpy trapz function Nadav Horesh wrote:
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
(Try not to top post on this list.) I can get it to work like this: import numpy as np y = np.arange(24).reshape(6,4) x = np.arange(6).reshape(-1,1) np.trapz(y, x, axis=0) From the text of the error message, you can see this is a problem with broadcasting. Due to broadcasting rules (which will *prepend* dimensions with size 1), you need to manually add an extra dimension to the end. Once I resize x, I can get this to work. You might want to look at this: http://www.scipy.org/EricsBroadcastingDoc 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
On Sat, Jul 12, 2008 at 10:30 PM, Nadav Horesh <nadavh@visionsense.com> wrote:
I am aware that the error is related to the broadcasting, and that it can be solved by matching the shape of x to that of y --- this is how I solved it in the first place. I was thinking that the function "promises" to integrate over an array given a x vector and the axis, so let obscure the broadcasting rules and just enable it to do the work. There is a reason to leave trapz as it is (or even drop it) since numpy should stay as close as possible to "bare metal", but this function is borrowed also by scipy integration package, thus I rather give it a face lift.
I think you can pretty much borrow the average function to do this, all you need to do is generate the proper weights and scaling. It's in numpy/lib/function_base.py Chuck
participants (3)
-
Charles R Harris -
Nadav Horesh -
Ryan May