function [pc,wpc,lambda,scores] = SpPC(t,x,q,wmu,dtyp)
%function [pc,wpc,lambda,scores] = SpPC(t,x,q,wmu,dtyp)
%
%Functional Spherical Principal Components
%
%INPUT: See comments on function SpMed.
% Input variables:
%   t: Time grid (as in SpMed).
%   x: Discretized curves (as in SpMed).
%   q: Number of PCs to estimate (scalar).
%   wmu: Spatial median weights. This is output 'w' of SpMed (n x 1 vector).
%   dtyp: Data type (as in SpMed; optional).
%
% Output variables:
%   pc: Spherical PCs (q x m matrix).
%   wpc: Weights (n x q matrix). (See paper for explanation).
%   lambda: Eigenvalues (q x 1 vector). (Note that these eigenvalues are
%   NOT consistent estimators of the true model eigenvalues; use e.g. median 
%   of squared PC scores to estimate eigenvalues consistently. See paper 
%   for further details).
%   scores: Individual PC scores (n x q matrix). Useful for computing
%   consistent estimators of the eigenvalues and for curve reconstruction
%   (x ~ ones(n,1)*med + scores*pc)
%
% External functions called: gsj.m
%

%% Input check
if nargin<4
    error('Not enough input variables')
elseif nargin==4
    dtyp = 's';
end
[n,m] = size(x);
if size(t,1)>1
    t = t';
end
if size(t,2)~=m
    error('Dimensions of T and X not compatible')
end
if size(wmu,2)>1
    wmu = wmu';
end
if size(wmu,1)~=n
    error('Dimensions of WMU and X not compatible')
end
if ~(strcmpi(dtyp,'n')||strcmpi(dtyp,'s'))
    error('Wrong input value for DTYP')
end

%% Inner-product matrices
A = 0.5*(x(:,1:m-1)*diag(t(2:m)-t(1:m-1))*x(:,1:m-1)' + ...
    x(:,2:m)*diag(t(2:m)-t(1:m-1))*x(:,2:m)');
if strcmpi(dtyp,'n')
    se2 = gsj(repmat(t',1,n),x');
    A1 = A - diag(se2)*(t(m)-t(1));
    if all(eig(A1)>-1e-6)
        A = A1;
    else
        disp('Corrected inner-product matrix is not nonnegative definite')
        disp('Using uncorrected matrix instead')
    end
end
B = (eye(n)-ones(n,1)*wmu')*A*(eye(n)-wmu*ones(1,n));
norms = sqrt(diag(B));
DN = zeros(n,1);
I = find(norms>eps);
DN(I) = 1./norms(I);
B = diag(DN)*B*diag(DN);
B = (B+B')/2;   % Just for numerical reasons

%% PC computation
opts.disp = 0;
[V,D] = eigs(B,q,'LA',opts);
D = diag(D);
if issorted(D(q:-1:1))
    wpc = V*diag(1./sqrt(D));
    lambda = D/n;
elseif issorted(D)
    wpc = V(:,q:-1:1)*diag(1./sqrt(D(q:-1:1)));
    lambda = D(q:-1:1)/n;
else
    [D,I] = sort(D,1,'descend');
    wpc = V(:,I)*diag(1./sqrt(D));
    lambda = D/n;
end
pc = wpc'*diag(DN)*(eye(n)-ones(n,1)*wmu')*x;
scores = 0.5*(pc(:,1:m-1)*diag(t(2:m)-t(1:m-1))*((eye(n)-ones(n,1)*wmu')*x(:,1:m-1))' + ...
    pc(:,2:m)*diag(t(2:m)-t(1:m-1))*((eye(n)-ones(n,1)*wmu')*x(:,2:m))');
scores = scores';