# HG changeset patch # User Rik # Date 1305392270 25200 # Node ID d8aff843a9e99cb93745ebb5742bd07d17874973 # Parent 68a59630798d1aae7e23003a7196232d897f58a8 Accept input x vector with y matrix for trapz,cumtrapz (bug #33292). * trapz.m, cumtrapz.m : Accept input vector, input matrix combination. Use diff() for 10% speedup. Add tests for ND-array operation. diff --git a/scripts/general/cumtrapz.m b/scripts/general/cumtrapz.m --- a/scripts/general/cumtrapz.m +++ b/scripts/general/cumtrapz.m @@ -50,53 +50,61 @@ print_usage (); endif - nd = ndims (x); - sz = size (x); + have_xy = have_dim = false; - have_x = false; - have_dim = false; if (nargin == 3) - have_x = true; + have_xy = true; have_dim = true; elseif (nargin == 2) if (! size_equal (x, y) && isscalar (y)) dim = y; have_dim = true; else - have_x = true; + have_xy = true; endif endif + if (have_xy) + nd = ndims (y); + sz = size (y); + else + nd = ndims (x); + sz = size (x); + endif + if (! have_dim) ## Find the first non-singleton dimension. - dim = find (sz > 1, 1); - if (isempty (dim)) - dim = 1; - endif + (dim = find (sz > 1, 1)) || (dim = 1); else - dim = floor (dim); - if (! (isscalar (dim) && 1 <= dim && dim <= nd)) - error ("cumtrapz: invalid dimension DIM"); + if (!(isscalar (dim) && dim == fix (dim)) + || !(1 <= dim && dim <= nd)) + error ("trapz: DIM must be an integer and a valid dimension"); endif endif n = sz(dim); - idx1 = cell (); - for i = 1:nd - idx1{i} = 1:sz(i); - endfor - idx2 = idx1; + idx1 = idx2 = repmat ({:}, [nd, 1]); idx1{dim} = 2 : n; idx2{dim} = 1 : (n - 1); - if (! have_x) + if (! have_xy) z = 0.5 * cumsum (x(idx1{:}) + x(idx2{:}), dim); else - if (! size_equal (x, y)) - error ("cumtrapz: X and Y must have the same shape"); + if (isvector (x) && !isvector (y)) + if (length (x) != sz(dim)) + error ("cumtrapz: length of X and length of Y along DIM must match"); + endif + ## Reshape vector to point along dimension DIM + shape = ones (nd, 1); + shape(dim) = sz(dim); + x = reshape (x, shape); + z = 0.5 * cumsum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim); + else + if (! size_equal (x, y)) + error ("cumtrapz: X and Y must have same shape"); + endif + z = 0.5 * cumsum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim); endif - z = 0.5 * cumsum ((x(idx1{:}) - x(idx2{:})) .* - (y(idx1{:}) + y(idx2{:})), dim); endif sz(dim) = 1; @@ -104,13 +112,23 @@ endfunction + %!shared x1,x2,y %! x1 = [0,0,0;2,2,2]; %! x2 = [0,2,4;0,2,4]; %! y = [1,2,3;4,5,6]; -%!assert (cumtrapz(y),[0,0,0;2.5,3.5,4.5]) -%!assert (cumtrapz(x1,y),[0,0,0;5,7,9]) -%!assert (cumtrapz(y,1),[0,0,0;2.5,3.5,4.5]) -%!assert (cumtrapz(x1,y,1),[0,0,0;5,7,9]) -%!assert (cumtrapz(y,2),[0,1.5,4;0,4.5,10]) -%!assert (cumtrapz(x2,y,2),[0,3,8;0,9,20]) +%!assert (cumtrapz(y), [0,0,0;2.5,3.5,4.5]) +%!assert (cumtrapz(x1,y), [0,0,0;5,7,9]) +%!assert (cumtrapz(y,1), [0,0,0;2.5,3.5,4.5]) +%!assert (cumtrapz(x1,y,1), [0,0,0;5,7,9]) +%!assert (cumtrapz(y,2), [0,1.5,4;0,4.5,10]) +%!assert (cumtrapz(x2,y,2), [0,3,8;0,9,20]) +%% Test ND-array implementation +%!shared x1,x2,y +%! x1 = 1:3; +%! x2 = reshape ([0,2,4;0,2,4], [1 2 3]); +%! y = reshape ([1,2,3;4,5,6], [1 2 3]); +%!assert (cumtrapz(y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3])) +%!assert (cumtrapz(x1,y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3])) +%!assert (cumtrapz(x2,y,3), reshape([0,3,8;0,9,20],[1 2 3])) + diff --git a/scripts/general/trapz.m b/scripts/general/trapz.m --- a/scripts/general/trapz.m +++ b/scripts/general/trapz.m @@ -64,60 +64,75 @@ print_usage (); endif - nd = ndims (x); - sz = size (x); + have_xy = have_dim = false; - have_x = false; - have_dim = false; if (nargin == 3) - have_x = true; + have_xy = true; have_dim = true; - endif - if (nargin == 2) + elseif (nargin == 2) if (! size_equal (x, y) && isscalar (y)) dim = y; have_dim = true; else - have_x = true; + have_xy = true; endif endif + if (have_xy) + nd = ndims (y); + sz = size (y); + else + nd = ndims (x); + sz = size (x); + endif + if (! have_dim) ## Find the first non-singleton dimension. - dim = find (sz > 1, 1); - if (isempty (dim)) - dim = 1; - endif + (dim = find (sz > 1, 1)) || (dim = 1); else - dim = floor (dim); - if (dim < 1 || dim > nd) - error ("trapz: invalid dimension DIM along which to sort"); + if (!(isscalar (dim) && dim == fix (dim)) + || !(1 <= dim && dim <= nd)) + error ("trapz: DIM must be an integer and a valid dimension"); endif endif n = sz(dim); - idx1 = cell (); - for i = 1:nd - idx1{i} = 1:sz(i); - endfor - idx2 = idx1; + idx1 = idx2 = repmat ({:}, [nd, 1]); idx1{dim} = 2 : n; idx2{dim} = 1 : (n - 1); - if (! have_x) + if (! have_xy) z = 0.5 * sum (x(idx1{:}) + x(idx2{:}), dim); else - if (! size_equal (x, y)) - error ("trapz: X and Y must have same shape"); + if (isvector (x) && !isvector (y)) + if (length (x) != sz(dim)) + error ("trapz: length of X and length of Y along DIM must match"); + endif + ## Reshape vector to point along dimension DIM + shape = ones (nd, 1); + shape(dim) = sz(dim); + x = reshape (x, shape); + z = 0.5 * sum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim); + else + if (! size_equal (x, y)) + error ("trapz: X and Y must have same shape"); + endif + z = 0.5 * sum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim); endif - z = 0.5 * sum ((x(idx1{:}) - x(idx2{:})) .* - (y(idx1{:}) + y(idx2{:})), dim); endif endfunction + %!assert (trapz(1:5), 12) %!assert (trapz(0:0.5:2,1:5), 6) -%!assert (trapz([1:5;1:5],2),[12;12]) -%!assert (trapz([1:5;1:5].',1),[12,12]) -%!assert (trapz([0:0.5:2;0:0.5:2],[1:5;1:5],2),[6;6]) -%!assert (trapz([0:0.5:2;0:0.5:2].',[1:5;1:5].',1),[6,6]) +%!assert (trapz([1:5;1:5].',1), [12,12]) +%!assert (trapz([1:5;1:5],2), [12;12]) +%!assert (trapz(repmat(reshape(1:5,1,1,5),2,2), 3), [12 12; 12 12]) +%!assert (trapz([0:0.5:2;0:0.5:2].',[1:5;1:5].',1), [6, 6]) +%!assert (trapz([0:0.5:2;0:0.5:2],[1:5;1:5],2), [6; 6]) +%!assert (trapz(repmat(reshape([0:0.5:2],1,1,5),2,2), ... +%! repmat(reshape(1:5,1,1,5),2,2), 3), [6 6; 6 6]) +%!assert (trapz(0:0.5:2,[(1:5)',(1:5)']), [6, 6]) +%!assert (trapz(0:0.5:2,[(1:5);(1:5)],2), [6; 6]) +%!assert (trapz(0:0.5:2,repmat(reshape(1:5,1,1,5),2,2),3), [6 6; 6 6]) +