function v = ppual(pp,x,left)
%PPUAL Evaluate function in ppform.
%
%   PPUAL(PP,X)  returns the value at X of the function whose ppform is in PP.
%
%   PPUAL(PP,X,'l<anything>')  takes the function to be left-continuous.
%   If the function in PP is m-variate and the third argument is actually an
%   m-cell, LEFT say, then, continuity from the left is enforced in the i-th 
%   variable if  LEFT{i}(1) = 'l'.  
%
%   The output is a matrix of size [d*m,n] if the function in PP is
%   univariate and d-vector valued, and [m,n] == size(X) .
%
%   If the function in PP is m-variate with m>1 and d-vector valued, then
%
%                         [d,n],         if X is of size [m,n]
%   the output is of size [d,n1,...,nm], if d>1  and X is {X1,...,Xm} 
%                         [n1,...,nm],   if d==1 and X is {X1,...,Xm} 
%
%   See also FNVAL, SPVAL, RSVAL, PPVAL.

%   Carl de Boor 19 jun 89
%   cb :  1 dec 1990 (correct misprint in last line)
%   cb : 19 sep 94 (temporarily programmed around the
%                   breaks([1 1 1 1 ... ]) discontinuity )
%   cb : 19 oct 94 (optional argument for left continuity)
%   cb :  9 may 95 (use .' instead of ', use reshape)
%   cb :  9 dec 95 (correct error in nargin==3 case and in handling matrix x )
%   cb : 15sep96 (handle empty X, cleaner if-statement; remove the
%                 fix for the former breaks([1 1 1 1 ..]) discontinuity,
%                 avoid multiplication when replicating a vector,
%                 return a matrix of size [d*m,n] with [m,n] = size(X) )
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   cb : 18oct97 (handle multivariate pp )
%   cb : 07may98 (admit former version of ppform; standardize and correct help)
%   cb : 04aug98 (slightly simplify formulation of max )
%   cb : 03sep99 (slightly updated the help)
%   cb : 23jun00 (enforce claim that only a LEFT of the form 'l...' will cause
%                 left continuity, particularly for scattered data sites;
%                 streamline determination of correct break interval)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.13 $

if ~isstruct(pp), pp = fn2fm(pp); end

if iscell(pp.breaks)   % we are dealing with a multivariate spline

   [breaks,coefs,l,k,d] = ppbrk(pp); m = length(breaks);

   if nargin>2 % set up left appropriately
      if ~iscell(left)
         temp = left; left = cell(1,m); [left{:}] = deal(temp);
      end
   else
      left = cell(1,m);
   end

   if iscell(x)  % evaluation on a mesh

      if length(x)~=m, error(['X should specify a(n) ',num2str(m),...
                              '-dimensional grid.']), end

      v = coefs; sizev = [d,l.*k]; % size(coefs); 
      nsizev = zeros(1,m);
      for i=m:-1:1
         nsizev(i) = length(x{i}(:)); dd = prod(sizev(1:m));
         v = reshape(ppual1(...
              ppmak(breaks{i},reshape(v,dd*l(i),k(i)),dd), x{i}, left{i} ),...
                [sizev(1:m),nsizev(i)]);
         sizev(m+1) = nsizev(i);
         if m>1
            v = permute(v,[1,m+1,2:m]); sizev(2:m+1) = sizev([m+1,2:m]);
         end
      end
      if d>1
         v = reshape(v,[d,nsizev]);
      else
         v = reshape(v,nsizev);
      end
 
   else          % evaluation at scattered points

      % locate the scattered data in the break sequences:
      [mx,n] = size(x);
      if mx~=m, error(['Each X(:,j) must be a ',num2str(m),'-vector.']), end
       
      ix = zeros(m,n);
      for i=1:m
         [ox,iindex] = sort(x(i,:));
         ix(i,iindex) = get_index(breaks{i},ox,left{i});
      end
      
      % ... and now set up lockstep polynomial evaluation
      % %%First,  select the relevant portion of the coefficients array. 
      % This has the additional pain that now there are k(i) coefficients
      % for the i-th univariate interval.  
      % The coefficients sit in the (m+1)-dimensional array COEFS, with
      % the (i+1)st dimension containing the coefficients in the i-th
      % dimension, and organized to have first the highest coefficients
      % for each interval, then the next-highest, etc (i.e., as if coming
      % from an array of size [l(i),k(i)]).
      % ix(:,j) is the index vector for the lower corner of j-th point
      % The goal is to extract, for the j-th point, the requisite coefficients
      % from the equivalent one-dimensional array for COEFS, computing a
      % base index from ix(:,j), and adding to this the same set of offsets
      % computed from the l(i) and k(i).
      
      temp = l(1)*[0:k(1)-1]'; 
      for i=2:m
         lt = length(temp(:,1));
         temp = [repmat(temp,k(i),1), ...
                 reshape(repmat(l(i)*[0:k(i)-1],lt,1),k(i)*lt,1)];
      end
      % also take care of the possibility that the function in PP is 
      % vector-valued 
      lt = length(temp(:,1));
      temp = [reshape(repmat([0:d-1].',1,lt),d*lt,1) temp(repmat(1:lt,d,1),:)];

      temp = num2cell(1+temp,1);
      offset = repmat(reshape(sub2ind(size(coefs),temp{:}),d*prod(k),1),1,n);
      
      temp = num2cell([ones(n,1) ix.'],1);
      base = repmat(sub2ind(size(coefs),temp{:}).',d*prod(k),1)-1;
      v = reshape(coefs(base+offset),[d,k,n]);
      
      % ... then do a version of local polynomial evaluation
      for i=m:-1:1
         s = reshape(x(i,:) - breaks{i}(ix(i,:)),[1,1,n]);
         otherk = d*prod(k(1:i-1));
         v = reshape(v,[otherk,k(i),n]);
         for j=2:k(i)
            v(:,1,:) = v(:,1,:).*s(ones(otherk,1),1,:)+v(:,j,:);
         end
         v(:,2:k(i),:) = [];
      end
      v = reshape(v,d,n);
   end
else
   if nargin<3, left = []; end
   v = ppual1(pp,x,left);
end

function v = ppual1(pp,x,left)
%PPUAL1 Evaluate univariate function in ppform.

[mx,nx] = size(x); lx = mx*nx; xs = reshape(x,1,lx);
%  if necessary, sort XS
tosort = 0;
if any(diff(xs)<0)
   tosort = 1; [xs,ix] = sort(xs);
end

%  take apart PP
[breaks,c,l,k,d] = ppbrk(pp);
%  if there are no points to evaluate at, return empty matrix of appropriate
%  size.
if lx==0, v = zeros(d,0); return, end

% for each data site, compute its break interval
index = get_index(breaks,xs,left);

% now go to local coordinates ...
xs = xs-breaks(index);
if d>1 % ... replicate XS and INDEX in case PP is vector-valued ...
   xs = reshape(xs(ones(d,1),:),1,d*lx);
   index = d*index; temp = [-d:-1].';
   index = reshape(1+index(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
end
% ... and apply nested multiplication:
v = c(index,1).';
for i=2:k
   v = xs.*v + c(index,i).';
end

v = reshape(v,d,lx);
if tosort>0,  v(:,ix) = v; end
v = reshape(v,d*mx,nx);

function ind = get_index(mesh,sites,left)
%GET_INDEX appropriate mesh intervals for given ordered data sites

if isempty(left)|left(1)~='l'
   ind = max(sorted(mesh(1:end-1),sites),1);
else
   ind = fliplr(max(length(mesh)-sorted(-fliplr(mesh(2:end)),-sites),1));
end

