function [xfit,fhat,fsig,dhat,dsig,ledf,spargcv]=SSfamfit(data,vspar,nfit,ab)
%function Smoothing Spline FAMily FIT
%Inputs
% data=[x,y]
% vspar=row vector of smoothing parameters
% nfit=number of design points where fitted values are computed
% ab=[a,b]=vector of endpoints of interval where fitted values are computed
%Outputs (can use any of 1st 2,3,5,6 or 7)
%    xfit is an nfit x 1 column vector,
%          locations where fitted values are calculated.
%    fhat is a nfit x length(vspar) matrix of smooths
%    fsig is a nfit x length(vspar) matrix of s.d. of fhat
%    dhat is a nfit x length(vspar) matrix of 1st deriv smooths
%    dsig is a nfit x length(vspar) matrix of s.d. of dhat
%    ledf is a nfit x length(vspar) matrix of local effective degrees of freedom
%    spargcv is an length(vspar) x 2 matrix, [vspar, gcv]
% By J.S. Marron and J.T. Zhang,   May 21, 2001


% Bin the data (allows handling large data sets, and uneven designs)
xmin = min(data(:,1)) ;
xmax = max(data(:,1)) ;
[bindat,bincent] = gplbinr(data,[xmin,xmax,nfit]) ;
flagn0 = bindat(:,1) > 0 ;
         %  1 where bin is nonempty
x = bincent(flagn0) ;
y = bindat(flagn0,2) ./ bindat(flagn0,1) ;
w = bindat(flagn0,1) ;
nbin = sum(flagn0) ;
         %  number of nonemty bins

% Bin the squares (for variance estimation)
bindat2 = gplbinr([data(:,1),data(:,2).^2],[xmin,xmax,nfit]) ;
y2s = bindat2(flagn0,2) ;
         % sum (NOT average) of Y's over each bin

flagw2 = w < 2 ;
nflagw2 = sum(flagw2) ;
if nflagw2 > 0 ;
  [temp,l] = sort(data(:,1)) ;
  yraw = data(l,2) ;
  nraw = size(yraw,1) ;
  oasig2 = sum((yraw(2:nraw) - yraw(1:(nraw-1))).^2) / (nraw - 1) ;
end ;


W=diag(w);
% diagonal elements are number of ties

%t=(x-xmin)/(xmax-xmin); % transform the data into [0,1] interval
t=x; %% no standardization

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating matrices Q,R, and K 
   h=t(2:nbin)-t(1:nbin-1);
   hinv=h.^(-1);

% Calculate Q:(nbin-1)x(nbin-2)
   Q=diag([hinv(1:nbin-2);0],-1)-diag([0;hinv(1:nbin-2)+hinv(2:nbin-1);0],0);
   Q=Q+diag([0;hinv(2:nbin-1)],1);
   Q=Q(2:nbin-1,:)';
   
% Calculate R:(nbin-2)x(nbin-2)
   R=diag(h(2:nbin-2),-1)/6+diag(h(1:nbin-2)+h(2:nbin-1))/3+diag(h(2:nbin-2),1)/6;

% Calculate K
   Rinv=inv(R);RinvQ=Rinv*Q';
   K=Q*RinvQ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%% Calculating nfit fits at ab=[a,b] which covers [t1,tn]
if nargin<4,
     a=.95*xmin;b=1.05*xmax; %% .05 can be selected
else
     a=ab(1);b=ab(2); %% ab=[a,b]=[left,right]
end          
xfit=a+(b-a)*([0:(nfit-1)]'/(nfit-1));%% nfit equ-spaced locations


%% Calculating coefficient matrices C,D and M for fitted values
%% at nfit equispaced points 
C=zeros(nfit,nbin);D=zeros(nfit,nbin-2);     
for i=1:nfit
    if xfit(i)<t(1),
         delta=(xfit(i)-t(1))*hinv(1);
         C(i,2)=delta;C(i,1)=1-delta;D(i,1)=delta*h(1)^2/6;
    elseif xfit(i)>=t(nbin),
         delta=(xfit(i)-t(nbin-1))*hinv(nbin-1);
         C(i,nbin)=delta;C(i,nbin-1)=1-delta;D(i,nbin-2)=C(i,nbin-1)*h(nbin-1)^2/6;
    elseif (xfit(i)>=t(1))&(xfit(i)<t(2)),
         delta=(xfit(i)-t(1))*hinv(1);
         C(i,2)=delta;C(i,1)=1-delta;
         D(i,1)=delta*(1-delta^2)*h(1)^2/6;     
    elseif (xfit(i)>=t(nbin-1))&(xfit(i)<t(nbin)),
         delta=(xfit(i)-t(nbin-1))*hinv(nbin-1);
         C(i,nbin)=delta;C(i,nbin-1)=1-delta;
         D(i,nbin-2)=delta*(1-delta)*(2-delta)*h(nbin-1)^2/6;
    else
       for j=2:(nbin-2),
          if (xfit(i)>=t(j))&(xfit(i)<t(j+1)),
                delta=(xfit(i)-t(j))*hinv(j);
                C(i,j+1)=delta;C(i,j)=1-delta;
                D(i,j)=delta*(1-delta^2)*h(j)^2/6;
                D(i,j-1)=delta*(1-delta)*(2-delta)*h(j)^2/6;
          end
       end         
    end
end               
M=C-D*RinvQ; %%%yfit=M*A*y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout >= 4 ;
  %% Calculating coefficient matrices C,D and M1 for first derivative
  %% fitted values at nfit  equispaced points 
  C=zeros(nfit,nbin);D=zeros(nfit,nbin-2);     
  for i=1:nfit
      if xfit(i)<t(1),
           C(i,2)=hinv(1);C(i,1)=-hinv(1);D(i,1)=h(1)/6;
      elseif xfit(i)>=t(nbin),
           C(i,nbin)=hinv(nbin-1);C(i,nbin-1)=-hinv(nbin-1);D(i,nbin-2)=-h(nbin-1)/6;
      elseif (xfit(i)>=t(1))&(xfit(i)<t(2)),
           delta=(xfit(i)-t(1))*hinv(1);
           C(i,2)=hinv(1);C(i,1)=-hinv(1);
           D(i,1)=(1-3*delta^2)*h(1)/6;
      elseif (xfit(i)>=t(nbin-1))&(xfit(i)<t(nbin)),
           delta=(xfit(i)-t(nbin-1))*hinv(nbin-1);
           C(i,nbin)=hinv(nbin-1);C(i,nbin-1)=-hinv(nbin-1);
           D(i,nbin-2)=(2-6*delta+3*delta^2)*h(nbin-1)/6;
      else
         for j=2:(nbin-2),
            if (xfit(i)>=t(j))&(xfit(i)<t(j+1)),
                  delta=(xfit(i)-t(j))*hinv(j);
                  C(i,j+1)=hinv(j);C(i,j)=-hinv(j);
                  D(i,j)=(1-3*delta^2)*h(j)/6;
                  D(i,j-1)=(2-6*delta+3*delta^2)*h(j)/6;
            end
         end         
      end
  end              
  M1=C-D*RinvQ; %%%dfit=M1*A*y
end ;

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
% Calculating fitted values at data points
%% selecting a sensible smoothing parameter by GCV

nspar=length(vspar);  
for i=1:nspar,
      %H=1+vspar(i)*diag(D);H=diag(H.^(-1));
      %A=Whinv*V*H*V'*Wh;
      A=pinv(W+vspar(i)*K)*W;
      yhat=A*y;
%% fits at design points

      %%% computing fhat and dhat
      MA=M*A;
      fhat(:,i)=MA*y;

      if nargout >= 3 ;

        error2=y2s - w.* y.^2;
        if nflagw2 > 0 ;  %  then replac local v. e. by global v.e.
          error2(flagw2) = w(flagw2) * oasig2 ;
        end ;
        sig2hat=A*error2;

        nhat=A*w;
        ysig2=sig2hat; %%ysig2=SS(error^2)/SS(w)
        fsig(:,i)=sqrt(diag(MA*diag(ysig2)*MA'));

        if nargout >= 4 ;
          M1A=M1*A;
          dhat(:,i)=M1A*y;
          dsig(:,i)=sqrt(diag(M1A*diag(ysig2)*M1A'));

         if nargout >= 6 ;
           oaedf = trace(A) ;
               %  overall effective degrees of freedom
           ledf(:,i) = oaedf ./ (M * nhat) ;
               %  local effective degrees of freedom

            %% GCV
            if nargout >= 7 ;
              gcv(i,:)=sum(w .* (y-yhat).^2)/(1-trace(A)/nbin)^2;
            end ;

          end ;

        end ;    

      end ;

 
end

if nargout >= 7 ;
  spargcv=[vspar',gcv];
end ;


