octave_packages/m/polynomial/ppval.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi})
  21 ## Evaluate the piecewise polynomial structure @var{pp} at the points @var{xi}.
  22 ## If @var{pp} describes a scalar polynomial function, the result is an
  23 ## array of the same shape as @var{xi}.
  24 ## Otherwise, the size of the result is @code{[pp.dim, length(@var{xi})]} if
  25 ## @var{xi} is a vector, or @code{[pp.dim, size(@var{xi})]} if it is a
  26 ## multi-dimensional array.
  27 ## @seealso{mkpp, unmkpp, spline, pchip}
  28 ## @end deftypefn
  29
  30 function yi = ppval (pp, xi)
  31
  32   if (nargin != 2)
  33     print_usage ();
  34   endif
  35   if (! (isstruct (pp) && strcmp (pp.form, "pp")))
  36     error ("ppval: first argument must be a pp-form structure");
  37   endif
  38
  39   ## Extract info.
  40   [x, P, n, k, d] = unmkpp (pp);
  41
  42   ## dimension checks
  43   sxi = size (xi);
  44   if (isvector (xi))
  45     xi = xi(:).';
  46   endif
  47
  48   nd = length (d);
  49
  50   ## Determine intervals.
  51   xn = numel (xi);
  52   idx = lookup (x, xi, "lr");
  53
  54   P = reshape (P, [d, n * k]);
  55   P = shiftdim (P, nd);
  56   P = reshape (P, [n, k, d]);
  57   Pidx = P(idx(:), :);#2d matrix size x: coefs*prod(d) y: prod(sxi)
  58
  59   if (isvector(xi))
  60     Pidx = reshape (Pidx, [xn, k, d]);
  61     Pidx = shiftdim (Pidx, 1);
  62     dimvec = [d, xn];
  63   else
  64     Pidx = reshape (Pidx, [sxi, k, d]);
  65     Pidx = shiftdim (Pidx, length (sxi));
  66     dimvec = [d, sxi];
  67   endif
  68   ndv = length (dimvec);
  69
  70   ## Offsets.
  71   dx = (xi - x(idx));
  72   dx = repmat (dx, [prod(d), 1]);
  73   dx = reshape (dx, dimvec);
  74   dx = shiftdim (dx, ndv - 1);
  75
  76   ## Use Horner scheme.
  77   yi = Pidx;
  78   if (k > 1)
  79     yi = shiftdim (reshape (Pidx(1,:), dimvec), ndv - 1);
  80   endif
  81
  82   for i = 2 : k;
  83     yi .*= dx;
  84     yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1);
  85   endfor
  86
  87   ## Adjust shape.
  88   if ((numel (xi) > 1) || (length (d) == 1))
  89     yi = reshape (shiftdim (yi, 1), dimvec);
  90   endif
  91
  92   if (isvector (xi) && (d == 1))
  93     yi = reshape (yi, sxi);
  94   elseif (isfield (pp, "orient") && strcmp (pp.orient, "first"))
  95     yi = shiftdim(yi, nd);
  96   endif
  97
  98   ##
  99   #if (d == 1)
 100   #  yi = reshape (yi, sxi);
 101   #endif
 102
 103 endfunction
 104
 105 %!shared b,c,pp,pp2,xi,abserr
 106 %! b = 1:3; c = ones(2); pp=mkpp(b,c);abserr = 1e-14;pp2=mkpp(b,[c;c],2);
 107 %! xi = [1.1 1.3 1.9 2.1];
 108 %!assert (ppval(pp,1.1), 1.1, abserr);
 109 %!assert (ppval(pp,2.1), 1.1, abserr);
 110 %!assert (ppval(pp,xi), [1.1 1.3 1.9 1.1], abserr);
 111 %!assert (ppval(pp,xi.'), [1.1 1.3 1.9 1.1].', abserr);
 112 %!assert (ppval(pp2,1.1), [1.1;1.1], abserr);
 113 %!assert (ppval(pp2,2.1), [1.1;1.1], abserr);
 114 %!assert (ppval(pp2,xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
 115 %!assert (ppval(pp2,xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
 116 %!assert (size(ppval(pp2,[xi;xi])), [2 2 4]);