Mercurial > hg > octave-nkf

changeset 6374:419017274c1e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-03-01 15:57:50 by jwe]
author jwe
date Thu, 01 Mar 2007 15:57:50 +0000
parents b3c37bc17c5f
children 5fced7a5eee8
files scripts/ChangeLog scripts/general/interp1.m
diffstat 2 files changed, 161 insertions(+), 46 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,7 @@
+2007-03-01  Paul Kienzle  <pkienzle@users.sf.net>
+
+	* general/interp1.m: Fix *style cases for decreasing x.
+
 2007-03-01  Muthiah Annamalai  <muthuspost@gmail.com>
 
 	* polynomial/roots.m: Check nargin before accessing arg.
--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -188,9 +188,9 @@
     endif
   elseif (strcmp (method, "*nearest"))
     if (pp)
-      yi = mkpp ([minx; minx+[0.5:(ny-1)]'*dx; maxx], y, szy(2:end));
+      yi = mkpp ([x(1); x(1)+[0.5:(ny-1)]'*dx; x(nx)], y, szy(2:end));
     else
-      idx = max (1, min (ny, floor((xi-minx)/dx+1.5)));
+      idx = max (1, min (ny, floor((xi-x(1))/dx+1.5)));
       yi(range,:) = y(idx,:);
     endif
   elseif (strcmp (method, "linear"))
@@ -210,10 +210,10 @@
   elseif (strcmp (method, "*linear"))
     if (pp)
       dy = [y(2:ny,:) - y(1:ny-1,:)];
-      yi = mkpp (minx + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
+      yi = mkpp (x(1) + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
     else
       ## find the interval containing the test point
-      t = (xi - minx)/dx + 1;
+      t = (xi - x(1))/dx + 1;
       idx = max(1,min(ny,floor(t)));
 
       ## use the endpoints of the interval to define a line
@@ -223,7 +223,7 @@
     endif
   elseif (strcmp (method, "pchip") || strcmp (method, "*pchip"))
     if (nx == 2 || method(1) == "*") 
-      x = linspace (minx, maxx, ny);
+      x = linspace (x(1), x(nx), ny);
     endif
     ## Note that pchip's arguments are transposed relative to interp1
     if (pp)
@@ -236,7 +236,7 @@
   elseif (strcmp (method, "cubic") || (strcmp (method, "*cubic") && pp))
     ## FIXME Is there a better way to treat pp return return and *cubic
     if (method(1) == "*")
-      x = linspace (minx, maxx, ny).'; 
+      x = linspace (x(1), x(nx), ny).'; 
       nx = ny;
     endif
 
@@ -277,7 +277,7 @@
 
     ## From: Miloje Makivic 
     ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
-    t = (xi - minx)/dx + 1;
+    t = (xi - x(1))/dx + 1;
     idx = max (min (floor (t), ny-2), 2);
     t = t - idx;
     t2 = t.*t;
@@ -293,7 +293,7 @@
 
   elseif (strcmp (method, "spline") || strcmp (method, "*spline"))
     if (nx == 2 || method(1) == "*") 
-      x = linspace(minx, maxx, ny); 
+      x = linspace(x(1), x(nx), ny); 
     endif
     ## Note that spline's arguments are transposed relative to interp1
     if (pp)
@@ -344,12 +344,27 @@
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
+## For each type of interpolated test, confirm that the interpolated
+## value at the knots match the values at the knots.  Points away
+## from the knots are requested, but only 'nearest' and 'linear'
+## confirm they are the correct values.
+
 %!shared xp, yp, xi, style
-%! xp=0:5;      yp = sin(2*pi*xp/5);
-%! xi = sort([-1, max(xp)*rand(1,6), max(xp)+1]);
+%! xp=0:2:10;      yp = sin(2*pi*xp/5);  
+%! xi = [-1, 0, 2.2, 4, 6.6, 10, 11];
+
 
+## The following BLOCK/ENDBLOCK section is repeated for each style
+##    nearest, linear, cubic, spline, pchip
+## The test for ppval of cubic has looser tolerance, but otherwise
+## the tests are identical.
+## Note that the block checks style and *style; if you add more tests
+## before to add them to both sections of each block.  One test, 
+## style vs. *style, occurs only in the first section.
+## There is an ENDBLOCKTEST after the final block
 %!test style = "nearest";
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1]), [NaN, NaN]);
+## BLOCK
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
@@ -358,11 +373,66 @@
 %!assert (isempty(interp1(xp,yp,[],style)));
 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 %!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
 %!assert (interp1(xp,[yp',yp'],xi,style),
-%!	  interp1(xp,[yp',yp'],xi,["*",style]));
-
-%!test style = "linear";
-%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
+%!test style=['*',style];
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
+## ENDBLOCK
+%!test style='linear';
+## BLOCK
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
+%!assert (interp1(xp,[yp',yp'],xi,style),
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
+%!test style=['*',style];
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
+## ENDBLOCK
+%!test style='cubic';
+## BLOCK
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
@@ -371,11 +441,32 @@
 %!assert (isempty(interp1(xp,yp,[],style)));
 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 %!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),100*eps);
+%!error interp1(1,1,1, style);
 %!assert (interp1(xp,[yp',yp'],xi,style),
 %!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
-
-%!test style = "cubic";
-%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
+%!test style=['*',style];
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),100*eps);
+%!error interp1(1,1,1, style);
+## ENDBLOCK
+%!test style='pchip';
+## BLOCK
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
@@ -384,11 +475,15 @@
 %!assert (isempty(interp1(xp,yp,[],style)));
 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 %!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
 %!assert (interp1(xp,[yp',yp'],xi,style),
-%!	  interp1(xp,[yp',yp'],xi,["*",style]),1000*eps);
-
-%!test style = "spline";
-%!assert (interp1(xp, yp, [-1, max(xp) + 1]), [NaN, NaN]);
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
+%!test style=['*',style];
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
 %!assert (interp1(xp,yp,xp,style), yp, 100*eps);
 %!assert (interp1(xp,yp,xp',style), yp', 100*eps);
 %!assert (interp1(xp',yp',xp',style), yp', 100*eps);
@@ -397,8 +492,47 @@
 %!assert (isempty(interp1(xp,yp,[],style)));
 %!assert (interp1(xp,[yp',yp'],xi(:),style),...
 %!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
+## ENDBLOCK
+%!test style='spline';
+## BLOCK
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
 %!assert (interp1(xp,[yp',yp'],xi,style),
-%!	  interp1(xp,[yp',yp'],xi,["*",style]),10*eps);
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
+%!test style=['*',style];
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,yp,xi,style),...
+%!	  interp1(fliplr(xp),fliplr(yp),xi,style),100*eps);
+%!assert (ppval(interp1(xp,yp,style,"pp"),xi),
+%!	  interp1(xp,yp,xi,style,"extrap"),10*eps);
+%!error interp1(1,1,1, style);
+## ENDBLOCK
+## ENDBLOCKTEST
 
 %!# test linear extrapolation
 %!assert (interp1([1:5],[3:2:11],[0,6],"linear","extrap"), [1, 13], eps);
@@ -407,7 +541,6 @@
 %!error interp1
 %!error interp1(1:2,1:2,1,"bogus")
 
-%!error interp1(1,1,1, "nearest");
 %!assert (interp1(1:2,1:2,1.4,"nearest"),1);
 %!error interp1(1,1,1, "linear");
 %!assert (interp1(1:2,1:2,1.4,"linear"),1.4);
@@ -424,25 +557,3 @@
 %!assert (interp1(1:2:8,1:2:8,1.4,"*cubic"),1.4);
 %!error interp1(1:2,1:2,1, "*spline");
 %!assert (interp1(1:2:6,1:2:6,1.4,"*spline"),1.4);
-
-%!assert (ppval(interp1(xp,yp,"nearest","pp"),xi),
-%!	  interp1(xp,yp,xi,"nearest","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"linear","pp"),xi),
-%!	  interp1(xp,yp,xi,"linear","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"cubic","pp"),xi),
-%!	  interp1(xp,yp,xi,"cubic","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"pchip","pp"),xi),
-%!	  interp1(xp,yp,xi,"pchip","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"spline","pp"),xi),
-%!	  interp1(xp,yp,xi,"spline","extrap"),10*eps);
-
-%!assert (ppval(interp1(xp,yp,"*nearest","pp"),xi),
-%!	  interp1(xp,yp,xi,"*nearest","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"*linear","pp"),xi),
-%!	  interp1(xp,yp,xi,"*linear","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"*cubic","pp"),xi),
-%!	  interp1(xp,yp,xi,"*cubic","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"*pchip","pp"),xi),
-%!	  interp1(xp,yp,xi,"*pchip","extrap"),10*eps);
-%!assert (ppval(interp1(xp,yp,"*spline","pp"),xi),
-%!	  interp1(xp,yp,xi,"*spline","extrap"),10*eps);