Mercurial > hg > octave-max

changeset 12667:d8aff843a9e9 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Rik <octave@nomad.inbox5.com>
date Sat, 14 May 2011 09:57:50 -0700
parents 68a59630798d
children e3dc23f7dd54
files scripts/general/cumtrapz.m scripts/general/trapz.m
diffstat 2 files changed, 91 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- 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]))
+
--- 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])
+