Re: Why is NaN used as extrapolation value?

[Søren Hauberg]

Søren Hauberg skrev:
> Anyway, I'll cook up a patch tonight unless somebody 
> gives me reason not to. Then we'll see what happens...
Okay, so I'm attaching a simple patch that replaces NaN with NA in the 
interpolation functions. I can't see how this could cause 
difficulties/incompatibilities, but as David points out, people do weird 
stuff. Apply or discard the patch as you will. The patch does fix one 
bug in 'interpn' where NaN's are being located by 'vi == NaN' (this 
should be 'isnan(vi)').

Søren

2007-06-18  Søren Hauberg   <address@hidden>
         * scripts/general/interp1.m scripts/general/interp2.m: replace
           NaN with NA.
         * scripts/general/interp3.m scripts/general/interpn.m: replace
           NaN with NA. Also change a comparison with NaN to a call to
           'isna'.
         * src/DLD-FUNCTIONS/__lin_interpn__.cc: replace octave_NaN with
           octave_NA.

Index: scripts/general/interp1.m
===================================================================
RCS file: /cvs/octave/scripts/general/interp1.m,v
retrieving revision 1.6
diff -u -r1.6 interp1.m
--- scripts/general/interp1.m   14 Jun 2007 06:56:42 -0000      1.6
+++ scripts/general/interp1.m   18 Jun 2007 11:09:14 -0000
@@ -51,7 +51,7 @@
 ##
 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
 ## the endpoints.  If @var{extrap} is a number, replace values beyond the
-## endpoints with that number.  If @var{extrap} is missing, assume NaN.
+## endpoints with that number.  If @var{extrap} is missing, assume NA.
 ##
 ## If the string argument 'pp' is specified, then @var{xi} should not be
 ## supplied and @code{interp1} returns the piece-wise polynomial that
@@ -95,7 +95,7 @@
   endif
 
   method = "linear";
-  extrap = NaN;
+  extrap = NA;
   xi = [];
   pp = false;
   firstnumeric = true;
@@ -378,7 +378,7 @@
 ## There is an ENDBLOCKTEST after the final block
 %!test style = "nearest";
 ## BLOCK
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -395,7 +395,7 @@
 %!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, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -412,7 +412,7 @@
 ## ENDBLOCK
 %!test style='linear';
 ## BLOCK
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -429,7 +429,7 @@
 %!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, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -446,7 +446,7 @@
 ## ENDBLOCK
 %!test style='cubic';
 ## BLOCK
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -463,7 +463,7 @@
 %!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, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -480,7 +480,7 @@
 ## ENDBLOCK
 %!test style='pchip';
 ## BLOCK
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -497,7 +497,7 @@
 %!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, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -514,7 +514,7 @@
 ## ENDBLOCK
 %!test style='spline';
 ## BLOCK
-%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NaN, NaN]);
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -531,7 +531,7 @@
 %!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, [min(xp)-1, max(xp)+1],style), [NA, NA]);
 %!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);
@@ -566,7 +566,7 @@
 %!error interp1(1,1,1, "*nearest");
 %!assert (interp1(1:2:4,1:2:4,1.4,"*nearest"),1);
 %!error interp1(1,1,1, "*linear");
-%!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"),[NaN,1,1.4,3,NaN]);
+%!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"),[NA,1,1.4,3,NA]);
 %!error interp1(1:3,1:3,1, "*cubic");
 %!assert (interp1(1:2:8,1:2:8,1.4,"*cubic"),1.4);
 %!error interp1(1:2,1:2,1, "*spline");
Index: scripts/general/interp2.m
===================================================================
RCS file: /cvs/octave/scripts/general/interp2.m,v
retrieving revision 1.7
diff -u -r1.7 interp2.m
--- scripts/general/interp2.m   12 Jun 2007 21:39:27 -0000      1.7
+++ scripts/general/interp2.m   18 Jun 2007 11:09:14 -0000
@@ -68,7 +68,7 @@
 ## If a scalar value @var{extrapval} is defined as the final value, then
 ## values outside the mesh as set to this value. Note that in this case 
 ## @var{method} must be defined as well. If @var{extrapval} is not
-## defined then NaN is assumed. 
+## defined then NA is assumed. 
 ##
 ## @seealso{interp1}
 ## @end deftypefn
@@ -92,7 +92,7 @@
 function ZI = interp2 (varargin)
   Z = X = Y = XI = YI = n = [];
   method = "linear";
-  extrapval = NaN;
+  extrapval = NA;
 
   switch (nargin)
     case 1
@@ -218,7 +218,7 @@
       ZI = Z(idx);
     endif
 
-    ## set points outside the table to NaN
+    ## set points outside the table to 'extrapval'
     ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval;
     ZI = reshape (ZI, shape);
   else
@@ -337,7 +337,7 @@
 %!  Orig = X.^2 + Y.^3;
 %!  xi = [0,4];
 %!  yi = [3,8]';
-%!  assert(interp2(x,y,Orig, xi, yi),[nan,nan;nan,nan]);
+%!  assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]);
 %!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);
 
 %!test % for values at boundaries
Index: scripts/general/interp3.m
===================================================================
RCS file: /cvs/octave/scripts/general/interp3.m,v
retrieving revision 1.2
diff -u -r1.2 interp3.m
--- scripts/general/interp3.m   14 Jun 2007 10:13:16 -0000      1.2
+++ scripts/general/interp3.m   18 Jun 2007 11:09:14 -0000
@@ -58,13 +58,13 @@
 ##
 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
 ## the endpoints.  If @var{extrap} is a number, replace values beyond the
-## endpoints with that number.  If @var{extrap} is missing, assume NaN.
+## endpoints with that number.  If @var{extrap} is missing, assume NA.
 ## @seealso{interp1, interp2, spline, meshgrid}
 ## @end deftypefn
 
 function vi = interp3 (varargin)
   method = "linear";
-  extrapval = NaN;
+  extrapval = NA;
   nargs = nargin;
 
   if (nargin < 1)
Index: scripts/general/interpn.m
===================================================================
RCS file: /cvs/octave/scripts/general/interpn.m,v
retrieving revision 1.3
diff -u -r1.3 interpn.m
--- scripts/general/interpn.m   14 Jun 2007 10:13:15 -0000      1.3
+++ scripts/general/interpn.m   18 Jun 2007 11:09:14 -0000
@@ -58,14 +58,14 @@
 ##
 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
 ## the endpoints.  If @var{extrap} is a number, replace values beyond the
-## endpoints with that number.  If @var{extrap} is missing, assume NaN.
+## endpoints with that number.  If @var{extrap} is missing, assume NA.
 ## @seealso{interp1, interp2, spline, ndgrid}
 ## @end deftypefn
 
 function vi = interpn (varargin)
 
   method = "linear";
-  extrapval = NaN;
+  extrapval = NA;
   nargs = nargin;
 
   if (nargin < 1)
@@ -150,7 +150,7 @@
   method = tolower (method);
   if (strcmp (method, "linear"))
     vi = __lin_interpn__ (x{:}, v, y{:});
-    vi (vi == NaN) = extrapval;
+    vi (isna(vi)) = extrapval;
   elseif (strcmp (method, "nearest"))
     yshape = size (y{1});
     yidx = cell (1, nd);
Index: src/DLD-FUNCTIONS/__lin_interpn__.cc
===================================================================
RCS file: /cvs/octave/src/DLD-FUNCTIONS/__lin_interpn__.cc,v
retrieving revision 1.1
diff -u -r1.1 __lin_interpn__.cc
--- src/DLD-FUNCTIONS/__lin_interpn__.cc        12 Jun 2007 21:39:27 -0000      
1.1
+++ src/DLD-FUNCTIONS/__lin_interpn__.cc        18 Jun 2007 11:09:14 -0000
@@ -188,7 +188,7 @@
   double *vi = Vi.fortran_vec ();
   octave_idx_type Ni = Vi.numel ();
 
-  double extrapval = octave_NaN;
+  double extrapval = octave_NA;
 
   for (int i = 0; i < n; i++)
     {