function y = polyvaln(p,x) 
%POLYVALN  evaluate a multivariate polynomial
%
%    The n-variate polynomial p, represented as an n-dimensional
%    array, is evaluated for each value in the array x. The last 
%    dimension of x should be of size n corresponding to the n
%    variables in the polynomial.
% 
%    Following the convention used in polyval(), the highest 
%    order coefficient of the polynomial comes first in the
%    array, and we take the first index to be the index of
%    the x-coefficient, the second of the y-coefficient and
%    so on. For a bivariate polynomial the array is a matrix 
%    and we have this sort of structure
%
%      [    x^n y^m,     x^n y^(m-1), ...     x^n y,    x^n
%       x^(n-1) y^m, x^(n-1) y^(m-1), ... x^(n-1) y, x^(n-1)
%                 :
%             x y^m,       x y^(m-1), ...        xy,      x
%               y^m,         y^(m-1), ...         y,      1]
%
%    Example: the polynomial x^2 + 2y^2 - 1 
%
%      p = [0,0,1;
%           0,0,0;
%           2,0,-1]
%
%    Evaluated at x=1, y=2 
%
%      polyvaln(p,[1,2])
%      ans = 8
% 
%    Evaluated at x,y = 1,..,5: we generate a 5x5x2 array xy
%    with the x values in xy(:,;,1), the y values in xy(:,:,2)
%
%      [x,y] = meshgrid(1:5)
%      xy = cat(3,x,y)
%      polyvaln(p,xy)
%      ans =
%          2    5   10   17   26
%          8   11   16   23   32
%         18   21   26   33   42
%         32   35   40   47   56
%         50   53   58   65   74
%
%    The evaluation method is a nested Horner's rule which is
%    implemented recursively. No attempt is made to exploit the 
%    sparsity (if any) of the polynomial. 
%
%    NOTE - this file should not be confused with the file of the
%    same name included with the PolyFit package by John D'Errico
%    http://mathworks.com/matlabcentral/fileexchange/34765-polyfitn
%
%    Copyright (c) 2011, 2012, J.J. Green
%    $Id: polyvaln.m,v 1.6 2012/05/15 20:26:02 jjg Exp $
%
%    Changes
%    
%    06 Jan 2011 : initial octave version
%    24 Jul 2011 : ported to work with matlab (7.12.0.635)
%    15 Mar 2012 : added note on the file in polyfit package

  % check arguments

  if nargin ~= 2
    error('exactly 2 arguments required')
  end

  % handle trivial case

  if isvector(p)
    y = polyval(p,x);
    return;
  end

  % size of last dimension of x

  if isvector(x)
    ndx  = 1;
    sldx = length(x);
  else
    dx   = size(x);
    ndx  = numel(dx);
    sldx = dx(ndx);
  end

  % number of dimensions of p

  dp = size(p);
  ndp = numel(dp);

  % check dimensions match

  if ndp ~= sldx
    error('size of last dim of x should equal number of dims of p')
  end

  % reshape x if needed

  if ndx > 1
    x = reshape(x,prod(dx(1:ndx-1)),sldx);
  end

  % call recursive subfunction

  y = pvn2(p,x);

  % reshape y to the same pattern as x if needed

  if ndx > 2

    y = reshape(y,dx(1:ndx-1));

  end

end

% recurse over the dimensions of p, evaluating using
% polyval() at the leaves of the recursion but pvn3()
% as we accumulate the totals

function y = pvn2(p,x) 

  if isvector(p)
    y = polyval(p,x);
    return
  end

  d  = size(p);
  nd = numel(d);

  c(1:nd-1) = {':'};
  x0 = x(:,2:nd);

  for i = 1:d(1)
    p0 = squeeze(p(i,c{:}));
    p1(:,i) = pvn2(p0,x0);
  end

  y = pvn3(p1,x(:,1));

end

% this is like polyval() except that p is an matrix
% with the same number of rows as x 

function y = pvn3(p,x) 

  n = size(p,2);

  y = p(:,1);
  for k = 2:n
    y = y.*x + p(:,k);
  end 

end
