function outstr = fksmean(t,x,order,ngrid,nkn,kn0)

%function outstr = fksmean(t,x,order,ngrid,nkn,kn0)
%Free-knot spline estimates of the mean of a sample of curves
%
%INPUT:
%   t       (1 x m)         Input time grid
%   x       (n x m)         Data matrix (subjects in rows)
%   order   (scalar)        Spline order
%   ngrid   (scalar)        Grid length
%   nkn     (1 x 2)         Min and max number of trial knots
%   kn0     (1 x nkn(1)-1)  Initial knot vector, if nkn(1)>1
%
%OUTPUT: Struct with 8 fields:
% These correspond to the p=nkn(2)-nkn(1)+1 knot sequences found by
% forward knot addition:
%   mu      (p x m)         Estimated means
%   knots   (1 x p cell)    Knot vectors
%   ase     (1 x p)         Average squared errors
%   gcv     (1 x p)         Generalized cross-validation criterion
% The fields mub, knotsb, aseb and gcvb are analogous and correspond to 
% knot sequences obtained by backward knot elimination, starting from 
% the optimum found by knot addition
%
%External calls: BSPL (B-spline basis), JUPP (Jupp transform of knots)

if nargin<6
    error('Not enough input arguments')
end
[n,m] = size(x);
if size(t,2)~=m
    error('Dimensions of X and T do not match')
end
p = nkn(2)-nkn(1)+1;
knots = cell(1,p);
grid = linspace(t(1),t(m),ngrid);
mu = zeros(p,m);
ase = Inf*ones(1,p);
gcv = Inf*ones(1,p);
options = optimset('Display','off','MaxIter',30,'LargeScale','off');

% Finding best knots for dimensions nkn(1) to nkn(2)
% Forward knot addition

for i = 1:p

    % Adds best knot from grid
    disp(' ')
    disp(['Number of knots = ' num2str(nkn(1)+i-1)])
    disp('Evaluating grid points ...')
    for j = 2:ngrid-1
        if i==1 & nkn(1)==1
            knots0 = grid(j);
        elseif i==1 & nkn(1)>1
            knots0 = sort([kn0 grid(j)]);
        else
            knots0 = sort([knots{i-1} grid(j)]);
        end
        B = bspl(t,order,[t(1) knots0 t(m)],0);
        mu0 = mean(x)*B*inv(B'*B)*B';
        ase0 = mean(mean((x-repmat(mu0,n,1)).^2));
        if ase0<ase(i)
            knots{i} = knots0;
            ase(i) = ase0;
        end
    end
    disp(['Best knots from grid search --> [' num2str(knots{i}) ']'])
    disp(['ASE --> ' num2str(ase(i))]);

    % Improves knots by local minimization
    jknots = jupp(knots{i},t(1),t(m),'direct');
    jknots = fminunc(@F,jknots,options,t,x,order);
    knots{i} = jupp(jknots,t(1),t(m),'inverse');

    % Evaluates ASE and GCV of improved knots
    B = bspl(t,order,[t(1) knots{i} t(m)],0);
    mu(i,:) = mean(x)*B*inv(B'*B)*B';
    ase(i) = mean(mean((x-repmat(mu(i,:),n,1)).^2));
    disp(['Adjusted knots --> [' num2str(knots{i}) ']'])
    disp(['ASE --> ' num2str(ase(i))])
    disp('--------------------')
    gcv(i) = ase(i)/(1-size(B,2)/(n*m))^2;
end

% Backward knot deletion

disp('**********')
disp(' ')
disp('Backward knot deletion')
[mcv,i0] = min(gcv);
knotsb = cell(1,i0-1);
mub = zeros(i0-1,m);
aseb = Inf*ones(1,i0-1);
gcvb = Inf*ones(1,i0-1);

for i = 1:i0-1
    
    % Finds least significant knot to delete
    ase0 = zeros(i0-i+1,1);
    knots0 = zeros(i0-i+1,i0-i);
    for j = 1:i0-i+1
        if i>1
            knots0(j,:) = knotsb{i-1}([1:j-1,j+1:i0-i+1]);
        else
            knots0(j,:) = knots{i0}([1:j-1,j+1:i0]);
        end
        B = bspl(t,order,[t(1) knots0(j,:) t(m)],0);
        mu0 = mean(x)*B*inv(B'*B)*B';
        ase0(j) = mean(mean((x-repmat(mu0,n,1)).^2));
    end
    if i>1
        [dif,j] = min(ase0-aseb(i-1));
    else
        [dif,j] = min(ase0-ase(i0));
    end
    aseb(i) = ase0(j);
    knotsb{i} = knots0(j,:);   
    disp(['Best knots after deletion --> [' num2str(knotsb{i}) ']'])
    disp(['ASE --> ' num2str(aseb(i))]);
    
    % Improves knots by local minimization
    jknots = jupp(knotsb{i},t(1),t(m),'direct');
    jknots = fminunc(@F,jknots,options,t,x,order);
    knotsb{i} = jupp(jknots,t(1),t(m),'inverse');
    
    % Evaluates ASE and GCV of improved knots
    B = bspl(t,order,[t(1) knotsb{i} t(m)],0);
    mub(i,:) = mean(x)*B*inv(B'*B)*B';
    aseb(i) = mean(mean((x-repmat(mub(i,:),n,1)).^2));
    disp(['Adjusted knots --> [' num2str(knotsb{i}) ']'])
    disp(['ASE --> ' num2str(aseb(i))])
    disp('--------------------')
    gcvb(i) = aseb(i)/(1-size(B,2)/(n*m))^2;
end

outstr = struct('mu',mu,'knots',{knots},'ase',ase,'gcv',gcv,...
    'mub',mub,'knotsb',{knotsb},'aseb',aseb,'gcvb',gcvb);


% ------------------

function ase = F(jknots,t,x,order)
[n,m] = size(x);
knots = jupp(jknots,t(1),t(m),'inverse');
B = bspl(t,order,[t(1) knots t(m)],0);
A = B'*B;
if rcond(A)<1e-15
    ase = mean(mean(x.^2));
else
    mu = mean(x)*B*inv(A)*B';
    ase = mean(mean((x-repmat(mu,n,1)).^2));
end