function [pc,lambda,scores] = RawPC(t,x,q,dtyp)
%function [pc,lambda,scores] = RawPC(t,x,q,dtyp)
%
%Functional Principal Components (no roughness penalization)
%
% Input variables:
%   t: Time grid (1 x m vector). Must be SORTED (in increasing order).
%   x: Discretized curves (n x m matrix). Each row represents a different
%      curve. Missing values NOT allowed. If your data set has missing
%      values, use either interpolation or smoothing to fill in the gaps.
%   q: Number of PCs to estimate (scalar).
%   dtyp: Data type (optional; character variable). Enter 's' if the data 
%         is "reasonably" smooth, or 'n' if the data is very noisy (in that
%         case a bias correction is necessary and the error variance is 
%         estimated nonparametrically using the GSJ estimator; this will
%         slow things down, so avoid if possible. It may also produce an 
%         inner-product matrix that is not nonnegative definite, in which
%         case the uncorrected matrix will be used). Default value: 's'.
%
% Output variables:
%   pc: Raw PCs (q x m matrix).
%   lambda: Eigenvalues (q x 1 vector).
%   scores: Individual PC scores (n x q matrix). Useful for individual
%       curve estimation (x ~ ones(n,1)*mean(x) + scores*pc)
%
% External functions called: gsj.m
%

%% Input check
if nargin<3
    error('Not enough input variables')
elseif nargin==3
    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 ~(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)/n)*A*(eye(n)-ones(n)/n);
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'*(eye(n)-ones(n)/n)*x;
scores = 0.5*(pc(:,1:m-1)*diag(t(2:m)-t(1:m-1))*((eye(n)-ones(n)/n)*x(:,1:m-1))' + ...
    pc(:,2:m)*diag(t(2:m)-t(1:m-1))*((eye(n)-ones(n)/n)*x(:,2:m))');
scores = scores';