function [out1,coefs,l,k,d] = ppbrk(pp,part)
%PPBRK Part(s) of a ppform.
%
%   [BREAKS,COEFS,L,K,D] = PPBRK(PP)  breaks the ppform in PP into its parts 
%   and returns as many of them as are specified by the output arguments. 
%
%   PPBRK(PP)  returns nothing, but prints all parts.
%
%   OUT1 = PPBRK(PP,PART)  returns the particular part specified by the string 
%   PART, which may be (the beginning character(s) of) one of the following
%   strings: 'breaks', 'coefs', 'pieces' or 'l', 'order' or 'k',
%   'dim'ension, 'interval'.
%   For a while, there is also the choice
%      'guide'
%   that returns the coefficient array in the form required in the ppform used
%   in `A Practical Guide to Splines', especially for PPVALU there. This is not
%   available for vector-valued and/or tensor product splines.
%
%   PJ = PPBRK(PP,J)  returns the ppform of the J-th polynomial piece of the 
%   function in PP.
%
%   PC = PPBRK(PP,[A B])  returns the restriction to the interval  [A .. B]  of 
%   the function in PP.
%
%   If PP contains an m-variate spline and PART is not a string, then it
%   must be a cell-array, of length m .
%
%   For example, if PP contains a bivariate spline with at least 4 pieces
%   in the first variable, then
%
%      ppp = ppbrk(pp,{4,[-1 1]});
%
%   gives the bivariate spline that agrees with the given one on the
%   rectangle  [pp.breaks{1}(4) .. [pp.breaks{1}(5)] x [-1 1] .
%
%   See also SPBRK, FNBRK.

%   Carl de Boor 25 feb 89
%   25 aug 94 (include request for just one item, incorporate pppce and ppcut,
%              simplify condition under which all the parts are printed.)
%   9 oct 95 (correct misprint), 19mar96 (clean up error(...))
%   20sep97 make certain it even works if PP is a column rather than a row.
%   cb : 18oct97 (change to structure, also cover multivariate spline)
%   cb : 08nov97 (also permit cut and piece for multivariate spline)
%   cb : 16nov97 (permit old array form, for backward compatibility)
%   cb : 05may98 (standardize the help)
%   cb : 09apr00 (provide helpful warning in case proposed basic interval
%                 is not proper, and avoid incomprehensible error return;
%                 also deal with possibility of trailing singleton dimension)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.11 $

if ~isstruct(pp)
   if pp(1)~=10
      error('The input array does not seem to describe a function in ppform.')
   else
      ppi = pp;
      di=ppi(2); li=ppi(3); bi=reshape(ppi(3+[1:li+1]),1,li+1); ki=ppi(5+li);
      pp.coefs=reshape(ppi(5+li+[1:di*li*ki]),di*li,ki);
      pp.form = 'pp'; pp.dim = di; pp.breaks = bi; pp.pieces = li; 
      pp.order = ki;
   end
end 

if (pp.form~='pp')
   error('The input does not seem to describe a function in ppform.')
end
if nargin>1
   if nargout>1, error('Too many output arguments for the given input.')
   end
   if isstr(part)
      switch part(1)
         case 'd',       out1 = pp.dim;
         case {'l','p'}, out1 = pp.pieces;
         case 'b',       out1 = pp.breaks;
         case {'k','o'}, out1 = pp.order;
         case 'c',       out1 = pp.coefs;
         case 'g',       % if the spline is univariate, nontensor product
                         % return the coefs in the form needed in the ppform
                         % used in PGS.
            if pp.dim>1|iscell(pp.order)
               error(sprintf(['''%s'' is only available for scalar-valued',...
                              ' univariate pp functions.'],part))
            else
               k = pp.order; cumprod([1 1:k-1].');
               out1 = (pp.coefs(:,k:-1:1).').*ans(:,ones(1,pp.pieces));
            end
         case 'i'
            if iscell(pp.breaks)
               for i=length(pp.order):-1:1
                  out1{i} = pp.breaks{i}([1 end]); end
            else
               out1 = pp.breaks([1 end]);
            end
         otherwise
            error(sprintf('''%s'' is not part of a ppform.',part))
      end
   else % we are to restrict PP to some interval or piece
      if iscell(part)  % we are dealing with a tensor-product spline

         [breaks,c,l,k,d] = ppbrk(pp); m = length(breaks);
         sizec = [d,l.*k]; %size(c);
         if length(sizec)~=m+1, error('Information in PP is inconsistent.'),
         end
         for i=m:-1:1
            dd = prod(sizec(1:m));
            ppi = ppbrk1(ppmak(breaks{i},reshape(c,dd*l(i),k(i)),dd),...
                        part{i}) ;
            breaks{i} = ppi.breaks; sizec(m+1) = ppi.pieces*k(i);
            c = reshape(ppi.coefs,sizec);
            if m>1
               c = permute(c,[1,m+1,2:m]);
               sizec(2:m+1) = sizec([m+1,2:m]);
            end
         end
         out1 = ppmak(breaks,c, sizec);

      else  % we are dealing with a univariate spline

         out1 = ppbrk1(pp,part);
      end
   end
else
   if nargout==0
     if iscell(pp.breaks) % we have a multivariate spline and, at present,
                          % I can't think of anything clever to do; so...
       disp(pp)
     else
       disp('breaks(1:l+1)'),        disp(pp.breaks)
       disp('coefficients(d*l,k)'),  disp(pp.coefs)
       disp('pieces number l'),      disp(pp.pieces)
       disp('order k'),              disp(pp.order)
       disp('dimension d of target'),disp(pp.dim)
     end
   else
      out1 = pp.breaks; coefs = pp.coefs; l = pp.pieces;
         k = pp.order;      d = pp.dim;
   end
end

function pppart = ppbrk1(pp,part)
%PPBRK1 restriction of pp to some piece or interval

if isempty(part), pppart = pp; return, end

[ignored, lp] = size(part);
if lp > 1 , % extract the part relevant to the interval 
            % specified by  part =: [a b]  
   pppart = ppcut(pp,part(1,1:2));
else        % extract the part-th polynomial piece of pp (if any)
   pppart = pppce(pp,part);
end

function ppcut = ppcut(pp,interv)
%PPCUT returns the part of pp  specified by the interval interv =: [a b]  

xl = interv(1); xr = interv(2); if xl>xr, xl = xr; xr = interv(1);  end
if xl==xr
   warning('No changes made since the given end points are equal.')
   ppcut = pp; return
end
 
%  the first pol. piece is  jl ,
% the one responsible for argument  xl
jl=pp.pieces; index=find(pp.breaks(2:jl)>xl); 
                                   % note that the resulting  index  ...
if (~isempty(index)), jl=index(1); % ... is shifted down by one  ...
end                                % ... because of  breaks(2: ...
%  if xl ~= breaks(jl), recenter the pol.coeffs.
x=xl-pp.breaks(jl);
di = pp.dim;
if x ~= 0
   a=pp.coefs(di*jl+[1-di:0],:);
   for ii=pp.order:-1:2
      for i=2:ii
         a(:,i)=x*a(:,i-1)+a(:,i);
      end
   end
   pp.coefs(di*jl+[1-di:0],:)=a;
end
 
%  the last pol. piece is  jr ,
% the one responsible for argument  xr .
jr=pp.pieces;index=find(pp.breaks(2:jr+1)>=xr); 
                                   % note that the resulting ...
if (~isempty(index)), jr=index(1); % index  is shifted down by
end                                % ... one because of  breaks(2: ...
 
%  put together the cut-down  pp
di = pp.dim;
ppcut = ppmak([xl pp.breaks(jl+1:jr) xr], ...
                        pp.coefs(di*(jl-1)+[1:di*(jr-jl+1)],:),di);

function pppce = pppce(pp,j)
%PPPCE returns the j-th polynomial piece of pp  (if any).

%  if  pp  has a  j-th  piece, ...
if (0<j)&(j<=pp.pieces)  %             ...  extract it
   di = pp.dim;
   pppce = ppmak([pp.breaks(j) pp.breaks(j+1)], ...
              pp.coefs(di*j+[1-di:0],:),di);
else
   error(['The given pp function does not have that many pieces.']);
end
