function spstr = spfpc(t,x,q,phi0,order,knots,lambda)
%function spstr = spfpc(t,x,q,phi0,order,knots,lambda)
%P-spline functional principal components 
%Computes q leading components orthogonal to phi0
%
%INPUT:
%   t       (1 x m)         Input time grid
%   x       (n x m)         Data matrix (subjects in rows)
%   q       (scalar)        Number of leading components to estimate
%   phi0    (q0 x m)        Initial principal components
%   order   (scalar)        Spline order
%   knots   (1 x p)         Knot vector (excluding interval endpoints)
%                           Repeated knots acceptable if lambda=0
%   lambda  (scalar)        Roughness penalty parameter
%
%OUTPUT: Struct with 3 fields:
%   pc      (q x m)         Estimated PCs (orthogonal to phi0)
%   l       (1 x q)         Eigenvalues associated with pc
%   scores  (n x q)         Individual component scores
%
%External calls: BSPL (B-spline basis)

if nargin<7
    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
if isempty(phi0)
    q0 = 0;
elseif size(phi0,2)~=m
    error('Dimensions of PHI0 and T do not match')
else
    q0 = size(phi0,1);
end
if order<=2 & lambda>0
    error('Order must be 3 or greater to penalize roughness')
end
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

B = bspl(t,order,[t(1) knots t(m)],0);
if lambda>0
    B2 = bspl(t,order,[t(1) knots t(m)],2);
else
    B2 = zeros(size(B));
end
[J,err] = chol(B'*diag(D)*B+lambda*B2'*diag(D)*B2);
if err>0 | any(isnan(J(:)))
    disp('Knots produce a degenerate spline regression matrix')
    spstr = [];
    return
end
p = length(knots)+order;
J1 = J\eye(p);

K = J1'*(B'*diag(D)*(x'*x/n)*diag(D)*B)*J1;
if q0>0
    Z = null(phi0*B*J1);
    K = Z'*K*Z;
end
[u1,d,u2] = svd(K);
if q0>0
    u1 = Z*u1;
end

pc = (B*J1*u1(:,1:q))';
norma = sqrt(pc*diag(D)*pc');
pc = pc/norma;
scores = x*diag(D)*pc';
l = mean(scores.^2);

spstr = struct('pc',pc,'l',l,'scores',scores);