function [g,gam,h,ssM,sM,dM] = smu_spline(x,y,alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Fit cubic smoothing spline for y as a function of x.
%
%  Use the Reinsch algorithm.
%  Algorithm based on the exposition in Green & Silverman (1994), sec 2.3
%      section 2.3
%
%  INPUTS:
%     x        n x 1 vector    Values of explanatory variable.
%                              Must be distinct and ordered x_i < x_{i+1}
%                              Require n >= 3
%     y        n x 1 vector    Values of response variable.
%     alpha    scalar          Value of smoothing parameter.
%                              Small for less smoothing, large for more smoothing.
%
%  OUTPUTS:
%     g         n x 1 vector   Values of spline at x:  g(x)
%     gam       (n-2) x 1 vector   Values of g"(x)
%                                  g"(x1) = 0 = g"(xn)
%     h         (n-1) x 1 vector  h_i = x_{i+1} - x_{i}
%     ssM       (n-2) x 1 vector  sub-subdiagonal of R + alpha*Q'Q      
%     sM        (n-2) x 1 vector  subdiagonal of R + alpha*Q'Q      
%     dM        (n-2) x 1 vector  diagonal of R + alpha*Q'Q      
%
%  Written by Karen Chiswell on 20 March 2006
%  Modified:
%       23 March 2006          Also output the bands of R + alpha*Q'Q
%        6 April 2006          Correct calculation of diagonals of R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = length(x);
nm1 = n-1;
nm2 = n-2;

%         initialize output
g = zeros(n,1);
gam = zeros(nm2,1);     %  first and last elements of gam are zero
h = zeros(nm1,1);

%----------------------------------------------------------------------------
%          Step 1:  compute h and Q'y
Qty = zeros(nm2,1);      %  note these elements are numbered 2:nm1
h(1) = x(2) - x(1);
for i = 2:nm1;
    h(i) = x(i+1) - x(i);
    Qty(i-1) = (y(i+1) - y(i)) / h(i) - (y(i) - y(i-1)) / h(i-1);
end   % loop on i

%----------------------------------------------------------------------------
%         Step 2:  compute bands of R + alpha*Q'Q and get cholesky factorization
%                  note that this is a (n-2)x(n-2) matrix with bandwidth 5
dQQ = zeros(nm2,1);       % diagonal of Q'Q
sQQ = zeros(nm2,1);       % subdiagonal of Q'Q  (last element is ignored)
ssQQ = zeros(nm2,1);      % sub-subdiagonal of Q'Q  (last 2 elements are ignored)
dR = zeros(nm2,1);        % diagonal of R
sR = zeros(nm2,1);        % subdiagonal of R  (last element is ignored)

for i=2:nm1;
    
    %    diagonal of R
    dR(i-1) = (h(i-1) + h(i))/3;
    %    subdiagonal of R
    if i < nm1
        sR(i-1) = h(i)/6;
    end
    
    %    diagonal of Q'Q
    dQQ(i-1) = (1/h(i-1))*(1/h(i-1)) ...
               + (-1/h(i-1)-1/h(i))*(-1/h(i-1)-1/h(i)) ...
               + (1/h(i))*(1/h(i));
    
    %    subdiagonal of Q'Q
    if i < nm1
        sQQ(i-1) = -(1/h(i-1) + 2/h(i) + 1/h(i+1)) / h(i);
    end
    
    %    sub-sub diagonal of Q'Q
    if i < (nm1-1)
        ssQQ(i-1) = (1/h(i+1)) * (1/h(i));
    end
    
end   % loop on i

%       put together pieces of R + alpha*Q'Q
dM = dR + alpha*dQQ;
sM = sR + alpha*sQQ;
ssM = alpha*ssQQ;

%       get Cholesky decomposition of the matrix R + alpha*Q'Q
[ssL,sL,d] = chol_band5(ssM,sM,dM);


%----------------------------------------------------------------------------
%         Step 3:  solve LDL' * gam = Q'y   for gam
gam = solve_band5(ssL,sL,d,Qty);


%----------------------------------------------------------------------------
%         Step 4:  g = y - alpha*Q*gam
Qg = zeros(n,1);
Qg(1) = gam(1)/h(1);
Qg(2) = gam(2)/h(2) - gam(1)*(1/h(1) + 1/h(2));
for i=3:nm2;
    Qg(i) = gam(i-2)/h(i-1) - gam(i-1)*(1/h(i-1) + 1/h(i)) + gam(i)/h(i);
end   % loop on i
Qg(nm1) = gam(nm2-1)/h(nm2) - gam(nm2)*(1/h(nm2) + 1/h(nm1));
Qg(n) = gam(nm2)/h(nm1);

g = y - alpha*Qg;

%    end of smu_spline.m