function outstr = fkspc(t,x,order,ngrid,nkn,kn0,phi0)
%function outstr = fkspc(t,x,order,ngrid,nkn,kn0,phi0)
%Free-knot spline estimates of functional principal components 
%Computes one PC at a time, orthogonal to phi0.
%
%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
%   phi0    (q x m)         Initial q principal components
%
%OUTPUT: Struct with 8 fields:
% These correspond to the p=nkn(2)-nkn(1)+1 knot sequences found by
% forward knot addition:
%   pc      (p x m)         Estimated leading PC (orthogonal to phi0)
%   knots   (1 x p cell)    Knot vectors
%   L       (1 x p)         Maximum eigenvalue associated with PC
%   cv      (1 x p)         Five-fold cross-validation criterion
% The fields pcb, knotsb, Lb and cvb are analogous and correspond to 
% knot sequences obtained by backward knot elimination, starting from 
% the optimum found by knot addition
%
%External calls: SPFPC (Spline functional PCs), JUPP (Jupp transform)

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

% Finding best knots for dimensions nkn(1) to nkn(2) by forward addition

for i = 1:p
    
    % Adds best knots from grid
    disp(' ')
    disp(['Number of internal 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
        pcstr = spfpc(t,x,1,phi0,order,knots0,0);
        if ~isempty(pcstr)
            if pcstr.l>L(i)
                knots{i} = knots0;
                L(i) = pcstr.l;
            end
        end
    end
    disp(['Best knots from grid search --> [' num2str(knots{i}) ']'])
    disp(['Max. eigenvalue --> ' num2str(L(i))]);
    
    % Improves knots by local maximization
    h = jupp(knots{i},t(1),t(m),'direct');
    h = fminunc(@F,h,options,t,x,order,phi0);
    knots{i} = jupp(h,t(1),t(m),'inverse');
    
    % Evaluates eigenvalue of PC for improved knots
    pcstr = spfpc(t,x,1,phi0,order,knots{i},0);
    if ~isempty(pcstr)
        phi(i,:) = pcstr.pc;
        L(i) = pcstr.l;
    end       
    disp(['Adjusted knots --> [' num2str(knots{i}) ']'])
    disp(['Max. eigenvalue --> ' num2str(L(i))])
    disp('--------------------')
end

% Five-fold cross-validation
if all(t(3:m)-t(2:m-1)==t(2)-t(1))
    D = ones(1,m)/m;
else
    D = [(t(1)+t(2))/2-t(1),(t(3:m)-t(1:m-2))/2,t(m)-(t(m-1)+t(m))/2];
end
cv = zeros(1,p);
ntest = ceil(n/5);
disp(' ')
for i = 1:p
    disp(['Cross-validating with ' num2str(nkn(1)+i-1) ' knots ...'])
    flag = 0;
    for j = 1:5
        itest = (j-1)*ntest+1:min(n,j*ntest);
        itrain = [1:(j-1)*ntest, min(n,j*ntest)+1:n];
        pcstr = spfpc(t,x(itrain,:),1,phi0,order,knots{i},0);
        if ~isempty(pcstr)
            s = x(itest,:)*diag(D)*pcstr.pc';
        else
            flag = 1;
            s = zeros(length(itest),1);
        end
        cv(i) = cv(i) + sum(s.^2); % We can ignore the PCs in phi0 here 
                                   % because they are orthogonal to new PC
    end
    cv(i) = cv(i)/n;
    if flag==1
        cv(i) = 0;
    end
end

% Backward knot deletion

disp('**********')
disp(' ')
disp('Backward knot deletion')
[mcv,i0] = max(cv);
knotsb = cell(1,i0-1);
phib = zeros(i0-1,m);
Lb = zeros(1,i0-1);

for i = 1:i0-1
    
    % Finds least significant knot to delete
    L0 = 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
        pcstr = spfpc(t,x,1,phi0,order,knots0(j,:),0);
        if ~isempty(pcstr)
            L0(j) = pcstr.l;
        end
    end
    if i>1
        [dif,j] = min(Lb(i-1)-L0);
    else
        [dif,j] = min(L(i0)-L0);
    end
    Lb(i) = L0(j);
    knotsb{i} = knots0(j,:);
    disp(['Best knots after deletion --> [' num2str(knotsb{i}) ']'])
    disp(['Max. eigenvalue --> ' num2str(Lb(i))]);

    % Improves knots by local maximization
    h = jupp(knotsb{i},t(1),t(m),'direct');
    h = fminunc(@F,h,options,t,x,order,phi0);
    knotsb{i} = jupp(h,t(1),t(m),'inverse');
    
    % Evaluates eigenvalue of PC for improved knots
    pcstr = spfpc(t,x,1,phi0,order,knotsb{i},0);
    if ~isempty(pcstr)
        phib(i,:) = pcstr.pc;
        Lb(i) = pcstr.l;
    end       
    disp(['Adjusted knots --> [' num2str(knotsb{i}) ']'])
    disp(['Max. eigenvalue --> ' num2str(Lb(i))])
    disp('--------------------')
end

% Five-fold cross-validation
ntest = ceil(n/5);
disp(' ')
cvb = zeros(1,i0-1);
for i = 1:i0-1
    disp(['Cross-validating with ' num2str(i0-i) ' knots ...'])
    cvb(i) = 0;
    flag = 0;
    for j = 1:5
        itest = (j-1)*ntest+1:min(n,j*ntest);
        itrain = [1:(j-1)*ntest, min(n,j*ntest)+1:n];
        pcstr = spfpc(t,x(itrain,:),1,phi0,order,knotsb{i},0);
        if ~isempty(pcstr)
            s = x(itest,:)*diag(D)*pcstr.pc';
        else
            flag = 1;
            s = zeros(length(itest),1);
        end
        cvb(i) = cvb(i) + sum(s.^2);
    end
    cvb(i) = cvb(i)/n;
    if flag==1
        cvb(i) = 0;
    end
end    

% Output
outstr = struct('pc',phi,'knots',{knots},'L',L,'cv',cv, ...
    'pcb',phib,'knotsb',{knotsb},'Lb',Lb,'cvb',cvb);

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

function y = F(h,t,x,order,phi0)
m = length(t);
knots = jupp(h,t(1),t(m),'inverse');
pcstr = spfpc(t,x,1,phi0,order,knots,0);
if ~isempty(pcstr), 
    y = -pcstr.l;
else
    y = 0; %Note that proper y's are always negative
end