Add interp1q

[David Bateman]

The benefits of the interp1q function are limited since with faster
machines than when this function was written for matlab, the portion of
time used to executed the error checking codein the interp1 function  is
typically very small. However, for completeness we should probably have
this function. Changeset attached that makes the same assumptions as the
matlab interp1q on the input parameters and just copies the relevant
code from interp1.m

D.

# HG changeset patch
# User David Bateman <address@hidden>
# Date 1210798279 -7200
# Node ID 60122ede5e7f1b6f91892d3ddbcef2a42164185e
# Parent  78adf97db4d8f73b45efa59271fa53340062b5d3
Add the interp1q function

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@ 2008-05-13  Bill Denney  <address@hidden
+2008-05-14  David Bateman  <address@hidden>
+
+       * general/interp1q.m: New function.
+       * general/Makefile.in (SOURCES): Add it here.
+
 2008-05-13  Bill Denney  <address@hidden>
 
        * general/isa.m: Use persistent cell arrays to hold class names
diff --git a/scripts/general/Makefile.in b/scripts/general/Makefile.in
--- a/scripts/general/Makefile.in
+++ b/scripts/general/Makefile.in
@@ -38,7 +38,7 @@ SOURCES = __isequal__.m __splinen__.m ac
   cart2sph.m cell2mat.m celldisp.m circshift.m common_size.m \
   cplxpair.m cumtrapz.m dblquad.m deal.m del2.m diff.m flipdim.m fliplr.m \
   flipud.m genvarname.m gradient.m idivide.m ind2sub.m int2str.m \
-  interp1.m interp2.m interp3.m interpn.m interpft.m \
+  interp1.m interp1q.m interp2.m interp3.m interpn.m interpft.m \
   is_duplicate_entry.m isa.m isdefinite.m isdir.m isequal.m \
   isequalwithequalnans.m isscalar.m issquare.m issymmetric.m \
   isvector.m logical.m logspace.m mod.m nargchk.m \
diff --git a/scripts/general/interp1q.m b/scripts/general/interp1q.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/interp1q.m
@@ -0,0 +1,70 @@
+## Copyright (C) 2008 David Bateman
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} address@hidden =} interp1q (@var{x}, @var{y}, 
@var{xi})
+## One-dimensional linear interpolation without error checking.
+## Interpolates @var{y}, defined at the points @var{x}, at the points
+## @var{xi}. The sample points @var{x} must be a strictly monotonically
+## increasing column vector. If @var{y} is an array, treat the columns
+## of @var{y} separately. If @var{y} is a vector, it must be a column
+## vector of the same length as @var{x}.
+##
+## Values of @var{xi} beyond the endpoints of the interpolation result
+## in NA being returned.
+##
+## Note that the error checking is only a significant portion of the
+## execution time of this @code{interp1} if the size of the input arguments
+## is relatively small. Therefore, the benefit of using @code{interp1q}
+## is relatively small.
+## @seealso{interp1}
+## @end deftypefn
+
+function yi = interp1q (x, y, xi)
+  x = x(:);
+  nx = size (x, 1);
+  szy = size (y);
+  ny = szy (1);
+  nc = prod (szy (2 : end));
+  y = reshape (y, ny, nc);
+  szx = size (xi);
+  xi = xi (:);
+  range = find (xi >= x (1) & xi <= x (nx));
+  yi = NA (size(xi, 1), size (y, 2));
+  xi = xi (range);
+  dy = y (2 : ny, :) - y (1 : ny - 1, :);
+  dx = x (2 : nx) - x (1 : nx - 1);
+  idx = lookup (x, xi, "lr");
+  s = (xi - x (idx)) ./ dx (idx);
+  yi (range, :) = s (:, ones (1, nc)) .* dy (idx, :) + y (idx, :);
+  if (length (szx) == 2 && any (szx == 1))
+    yi = reshape (yi, [max (szx), szy (2 : end)]);
+  else
+    yi = reshape (yi, [szx, szy(2:end)]);
+  endif
+endfunction
+
+%!shared xp, yp, xi
+%! xp=[0:2:10].';      yp = sin(2*pi*xp/5);
+%! xi = [-1; 0; 2.2; 4; 6.6; 10; 11];
+%!assert (interp1(xp, yp, [min(xp)-1; max(xp)+1]), [NA; NA]);
+%!assert (interp1(xp,yp,xp), yp, 100*eps);
+%!assert (isempty(interp1(xp,yp,[])));
+%!assert (interp1(xp,[yp,yp],xi), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);
+%!assert (interp1(xp,yp,[xi,xi]), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);
+%!assert (interp1(xp,[yp,yp],[xi,xi]), [interp1(xp,yp,xi),interp1(xp,yp,xi)]);