function pp = ppmak(breaks,coefs,d)
%PPMAK Put together a spline in ppform.
%
%   PPMAK(BREAKS,COEFS)  puts together a spline in ppform from the breaks
%   BREAKS and coefficient matrix COEFS. Each column of COEFS is
%   taken to be one coefficient, i.e., the spline is taken to be D-vectorvalued
%   if COEFS has D rows. Further, with L taken as  length(BREAKS)-1, the order
%   K of the spline is computed as (# cols(COEFS))/L, and COEFS is
%   interpreted as a three-dimensional array of size [D,K,L], with 
%   COEFS(:,:,j) containing the local polynomial coefficients for the j-th
%   polynomial piece, from highest to lowest.
%
%   PPMAK  will prompt for BREAKS and COEFS.
%
%   PPMAK(BREAKS,COEFS,D) expects the matrix COEFS to be of size [D*L,K] and
%   so that, taken as a three-dimensional array of size [D,L,K], it has the 
%   local polynomial coefficients of the j-th polynomial piece, from highest 
%   to lowest, in COEFS(:,j,:). In particular, the order K is taken to be the 
%   column number of COEFS, and L is taken to be (# rows(COEFS))/D, and
%   BREAKS is expected to be of length L+1.
%   The toolbox uses internally only this second format.
%
%   For example,  ppmak([1 3 4],[1 2 5 6;3 4 7 8])  and
%                 ppmak([1 3 4],[1 2;3 4;5 6;7 8],2)
%   specify the same function (2-vector-valued, of order 2).
%   The toolbox uses internally only the second format.
%
%   PPMAK({BREAKS1,...,BREAKSm},COEFS)  puts together an m-variate 
%   tensor-product spline in ppform. In this case, COEFS is expected to be of 
%   size [D,lk], with  lk := [l1*k1,...,lm*km]  and  li = length(BREAKS{i})-1, 
%   all i , and this defines  D  and  k := [k1,...,km].  If, instead, COEFS is 
%   only an m-dimensional array, then  D  is taken to be 1. 
%
%   PPMAK({BREAKS1,...,BREAKSm},COEFS,SIZEC)  uses the optional third argument
%   to specify the size of COEFS explicitly in case one or more of its trailing
%   dimensions is a singleton and thus appears to be of lower dimension.
%   
%   For example, if we intend to construct a 2-vector-valued bivariate 
%   polynomial on the rectangle [-1 .. 1] x [0 .. 1], linear in the first
%   variable and constant in the second, say
%      coefs = zeros(2,2,1); coefs(:,:,1) = [1 0; 0 1];
%   then the straightforward
%      pp = ppmak({[-1 1],[0 1]},coefs);
%   will fail, producing a scalar-valued function of order 2 in each variable,
%   as will 
%      pp = ppmak({[-1 1],[0 1]},coefs,size(coefs));
%   while the command
%      pp = ppmak({[-1 1],[0 1]},coefs,[2 2 1]);
%   will succeed.
%
%   See also PPBRK, RPMAK, SPMAK, RSMAK, FNBRK.

%   Carl de Boor 1 jun 89
%   cb : jun 94 (change from  error('')  to  error(' ') )
%   cb :  9 may '95 (use .' instead of '); 3 nov 95 (update help)
%   matlab: 19mar96 (clean up error(...))
%   cb : 15aug96 (test for nontrivial basic interval)
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   cb : 18oct97 (change to structure, also cover multivariate spline)
%   cb : 07dec97 (reorder form for backward consistency, ensure breaks is 1-row)
%   cb : 07may98 (update and standardize help)
%   cb : 23jan00 (skip warning about third argument in multivariate case)
%   cb : 09apr00 (deal with fact that matlab drops trailing singleton dims)
%   cb : 23jun00 (... and improve the help by adding an example re this)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.13 $

if nargin==0
   breaks=input('Give the (l+1)-vector of breaks  >');
   coefs=input('Give the (d by (k*l)) matrix of local pol. coefficients  >');
end

if iscell(breaks)  % we are dealing with a tensor-product spline
   m = length(breaks); sizec = size(coefs); 
   if nargin>2
      if prod(sizec)~=prod(d)
        error('The coefficient array is not of the explicitly specified size.')
      end, sizec = d;  end
   switch length(sizec)
      case m  % coefficients of a scalar-valued function
         sizec = [1 sizec]; coefs = reshape(coefs, sizec); 
      case m+1
      otherwise
         error(['If BREAKS is a cell-array of length m, then COEFS must ',...
                'be an (m)- or (m+1)-dimensional array.'])
   end
   d = sizec(1);
   for i=m:-1:1
      l(i) = length(breaks{i})-1;
      k(i) = fix(sizec(i+1)/l(i));
      if k(i)<=0|k(i)*l(i)~=sizec(i+1)
         error(sprintf(['The specified number ',num2str(l(i)),' of',...
         ' polynomial pieces is incompatible\nwith the total number ',...
         num2str(sizec(i+1)),' supplied in variable ',num2str(i),'.']))
      end
      breaks{i} = reshape(breaks{i},1,l(i)+1);
   end
else
  if nargin<3
     [d,kl]=size(coefs);
     if isempty(coefs)
        error('The coefficient sequence is empty.')
     end
     l=length(breaks)-1;k=fix(kl/l);
     if (k<=0)|(k*l~=kl)
       error(sprintf(['The specified number %.0f of polynomial pieces is',...
        ' incompatible\nwith the total number %.0f of coefficients',...
        ' supplied.'],l,kl));
        elseif any(diff(breaks)<0)
        error('The break sequence should be nondecreasing.')
     elseif breaks(1)==breaks(l+1)
        error('The extreme breaks should be different.')
     else
        % the ppformat expects coefs in array  (d*l) by k, while the standard
        % input supplies them in an array d by (k*l) . This requires the
        % following shuffling, from  D+d(-1+K + k(-1+L))=D-d +(K-k)d + dkL
        % to  D+d(-1+L + l(-1+K)=D-d +(L-l)d + dlK .
        % This used to be handled by the following:
        % c=coefs(:); temp = ([1-k:0].'*ones(1,l)+k*ones(k,1)*[1:l]).';
        % coefs=[1-d:0].'*ones(1,kl)+d*ones(d,1)*(temp(:).');
        % coefs(:)=c(coefs);
        % Thanks to multidimensional arrays, we can now simply say
        coefs = reshape(permute(reshape(coefs,[d,k,l]),[1,3,2]),d*l,k);
     end
  else
     [dl,k] = size(coefs); l = dl/d;
     if l+1~=length(breaks)
        error('Input arguments have incompatible sizes.'), end
  end
  breaks = reshape(breaks,1,l+1);
end
pp.form = 'pp';
pp.breaks = breaks;
pp.coefs = coefs;
pp.pieces = l;
pp.order = k;
pp.dim = d;
