function [xfit,fhat,fsig,dhat,dsig,ledf,spargcv,vedf]=SSfamfit0(data,vspar,nfit,ab,indfig)
%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
% indfig=1  draw pictures
%%      =0   no pictures
%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]
%    vedf  vector of effective degrees of freedom
% By J.S. Marron and J.T. Zhang,   May 21, 2001

if nargin<5|length(indfig)==0,
   indfig=0;
end   

if nargin<3|length(nfit)==0, 
    nfit=min([length(data(:,1)),401]);
    %nfit=301;
end    

%% sort the data
[x,temp]=sort(data(:,1));
data=data(temp,:);

% 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

%%% Data Transformation

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

if ~(nargin<2|length(vspar)==0)
   vspar=vspar*xrange^4;
end   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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]
%% [a,b] is tansferred into [0,1]
if nargin<4|length(ab)==0,
     a=0;b=1; %% .05 can be selected
else
     a=(ab(1)-xmin)/xrange;b=(ab(2)-xmin)/xrange; %% 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

Wh=diag(w.^(.5));
Whinv=diag(w.^(-.5));
G=Whinv*K*Whinv;
[U,D]=eig(G);
[d,temp]=sort(abs(diag(D)));
%U=U(:,temp);
%UL=Whinv*U;UR=U'*Wh;

if nargin<2 | length(vspar)==0, %%specify a range of vspar
   
  dmin=d(3);%% First 3 could be 0     
  dmax=d(nbin);
  sparmin=(1/.8-1)/dmax;
  sparmax=(nbin/2.2-1)/dmin;
  vspar=logspace(log10(sparmin),log10(sparmax),11); %% 11 spars
  
end

nspar=length(vspar);  
for i=1:nspar,

      %% fits at design points
      %A=UL*diag((1+vspar(i)*d).^(-1))*UR;
      A=pinv(W+vspar(i)*K)*W;
      if nargout>=8,
        vedf(i)=trace(A); %% added by Jin-Ting Zhang  Jan. 15, 2002.
      end
      yhat=A*y;
 
      %%% 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 ;
              %z=UR*y;
              %temp=nbin*vspar(i)*d./(1+vspar(i)*d);
              %gcv(i,:)=sum(temp.^2.*z.^2)/(mean(temp))^2;
              gcv(i,:)=sum(w.*(y-yhat).^2./(1-diag(A)).^2);
            end ;

          end ;

        end ;    

      end ;

 
end

if nargout >= 7 ;
  spargcv=[vspar(:)/xrange^4,gcv];
end ;


xfit=xmin+xrange*xfit;


if indfig==1&nargout>=2,
  figure(1),clf
  plot(data(:,1),data(:,2),'o')
  hold on
  plot(xfit,fhat)
  hold off
end  

if indfig==1& nargout>=4,
  figure(2),clf
  plot(xfit,dhat)
end  

if indfig==1 & nargout>=7,
   figure(3),clf
   plot(log10(spargcv(:,1)),spargcv(:,2),'-o')
 end