function [yhat,dhat,hgcv,vhgcv]=monofit(data,params)
%function [yhat,dhat,hgcv,vhgcv]=monofit(data,params)
%% Inputs  
%%    data=[x,y] nx2 matrix
%%    params=[h,mflag,method,indfig]
%%             h     smoothing parameter, 
%%           if not provided, selected by GCV
%%            REMARK: the h selected by GCV may occasionally be ridicules
%%       mflag = 0   regular fit with no constraint
%%             = 1  monotone increasing
%%             = 2  montone  decreasing
%%       method = 1  Method 1 
%%              = 2  Method 2
%%       indfig = 0 (default) no plots for fits
%%              = 1  plots of fits, 1st derivative fits, GCV, SSE
%% Outputs
%%    yhat=[t,yhat,ysig]
%%    dhat=[t,ghat,gsig] 
%%    hgcv=[h,gcv]
%     vhgcv=[vh,vgcv]
%  
% By J.T. Zhang, Jan.11,1998.
% Revised  April 29, 1999
%% Revised 13/12/2001, Jin-Ting Zhang at NUS
%% Revised Oct. 3, 2005

 
[temp,flag]=sort(data(:,1));
data=data(flag,:);   
x=data(:,1);y=data(:,2);n=size(x,1);
[ux,flag1,flag2]=unique(x);
xmin=min(ux);xmax=max(ux);nux=length(ux);

if nargin<2|length(params)==0,
  spar=0;mflag=0;method=2;indfig=0;
elseif length(params)==1,
  spar=params(1);mflag=0;method=2;indfig=0;
elseif length(params)==2,
  spar=params(1);mflag=params(2);method=2;indfig=0;
elseif length(params)==3,
  spar=params(1);mflag=params(2);method=params(3);indfig=0; 
elseif length(params)==4,
  spar=params(1);mflag=params(2);method=params(3);indfig=params(4);   
end  

%% Indicine matrix, size n x nux
E=((x*ones(1,nux))==(ones(n,1)*ux'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculating the K matrix: (nux-1) x (nux-1))
   nn=n;n=nux;
   h=ux(2:n)-ux(1:(n-1));
   hinv=h.^(-1);

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

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

%% Compute the monotone-smoothing related  matrices  
bK=[[0,zeros(1,n)];[zeros(n,1),K]]; % forming a big matrix

if method==1,

% Calculating the matrix M
%%S1: nxn
 S1(1,:)=zeros(1,n);
 for i=1:(n-1),
  S1(i+1,:)=[h(1:i)',zeros(1,n-i)]+[0,h(1:i)',zeros(1,n-i-1)];
 end  

%%S2: nx(n-2)
 S2(1,:)=zeros(1,n-2);
 S2(2,:)=[h(1),zeros(1,n-3)];
 for i=1:(n-2),
   S2(i+2,:)=[h(1:i)',zeros(1,n-i-2)]+[h(2:i)',zeros(1,n-i-1)];
 end  
  S2=S2.^3;

  M=S1/2-S2/24*Rinv*Q';
  bM=[ones(n,1),M];% form a big M: nx(n+1)
end

% Calculating the matrix S: n x (n+1)
S(1,:)=[1,zeros(1,n)];
for i=1:(n-1),
  S(i+1,:)=[1,h(1:i)'/2,zeros(1,n-i)]+[0,0,h(1:i)'/2,zeros(1,n-1-i)];
end 
n=nn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%  

if method==1,
   X=E*bM;X0=E*S;
else
   X=E*S;
end
H0=X'*X;b=X'*y;      

if spar==0,

  [U,D]=eig(H0);d=abs(diag(D)); 
  B=U*pinv(diag(d.^(.5)))*U';
  
  [U,D]=eig(B*bK*B');d=sort(abs(diag(D)));
  hmin=((nux+1)/(nux-6)-1)/max(d);
  hmax=((nux+1)/2.1-1)/d(3);
  
  for ii=1:3,
  
  format short e 
  nspar=5;
  vspar=logspace(log10(hmin),log10(hmax),nspar);
  
  for i=1:nspar,
    spar=vspar(i);
    H=pinv(H0+spar*bK);
    beta=H*b; %%
    yhat0=beta(1);
    ghat=beta(2:(nux+1)); %% first derivatives
    
    if mflag==0, %%no constraits
       yhat=X*beta;
    elseif mflag==1, %% monotone increasing 
       yhat=E*S*[yhat0;ghat.*(ghat>=0)];     
    elseif mflag==2, %% monotone decreasing
       yhat=E*S*[yhat0;ghat.*(ghat<=0)];
    end
    
    if method==1,
        df=trace(H*X'*X0);
    elseif method==2,
        df=trace(H*X'*X);
    end        
    sse=sum((y-yhat).^2);
    gcv=sse/(1-df/n)^2;
    vhgcv((ii-1)*nspar+i,:)=[spar,gcv,sse,df];
  end  

  [temp,flag0]=min(vhgcv(:,2));
  spar=vhgcv(flag0(1),1);
  hmin=spar*10^(-2);hmax=spar*10^3;
  
end  

[temp,flag]=sort(vhgcv(:,1));
vhgcv=vhgcv(flag,:);
  
  
else
  nspar=1;  
end  

%% fitting
    H=pinv(H0+spar*bK);
    beta=H*b; %%differences starting from 4th entry
    yhat0=beta(1);
    ghat=beta(2:(nux+1)); %% first derivatives
    
   
   if mflag==0, %%no constraits
       yhat=X*beta;
    elseif mflag==1, %% monotone increasing  
          yhat=E*S*[yhat0;ghat.*(ghat>=0)];
    elseif mflag==2, %% monotone decreasing  
          yhat=E*S*[yhat0;ghat.*(ghat<=0)];
    end    
  
   if method==1,
     df=trace(H*X'*X0);
   elseif method==2,  
     df=trace(H*X'*X);
   end  
  sse=sum((y-yhat).^2);
  gcv=sse/(1-df/n)^2;
  hsig2=sse/(n-df);
  
  if method==1,
     A=X0*H*X';
  elseif method==2,
     A=X*H*X';
  end
        
  ysig=sqrt(hsig2*diag(A*A'));
  
  B=H*X';B=B(2:(nux+1),:);
  gsig=sqrt(hsig2*diag(B*B'));
 
  hgcv=[spar,gcv,sse,df];
  
  if nspar==1, vhgcv=hgcv;end 
    
if indfig==1  % plotting
  subplot(2,1,1)
  % Drawing a plot
  plot(data(:,1),data(:,2),'o',x,yhat,'r-',...
      x,[yhat-1.96*ysig,yhat+1.96*ysig],'g--')
 
  xlabel('x'),ylabel('y'),title('(a) Monotone Function Fit')
 
  
  subplot(2,2,3)
  plot(ux,ghat,'r-',...
       ux,[ghat-1.96*gsig,ghat+1.96*gsig],'g--')
  title('(b) 1st Derivative Fit')
  
end

if nspar>1,
  
  subplot(2,2,4)
  plot(log10(vhgcv(:,1)),log10(vhgcv(:,2)),'-o')
  xlabel('log10(\lambda)')
  ylabel('log10(GCV)')
  title('(c) GCV Curve')
 
end


yhat=[x,yhat,ysig];
dhat=[ux,ghat,gsig];