function [r,w,mu,pc,l,sc] = fte(t,x,q,alpha,beta,weight)
%function [r,w,mu,pc,l,sc] = fte(t,x,q,alpha,beta,weight)
%Functional Trimmed Estimators (for univariate, vector-valued functions)
%
%INPUT:
%   T       (1 x m)         Time grid
%   X       (n x m x p)     Data curves (p dimensional, n replications, sampled at T)
%   Q       scalar          Number of principal components to be computed (Q>=0)
%   ALPHA   scalar          Radii are ALPHA-quantiles of the distances (usually .50)
%   BETA    scalar          Trimming proportion for weights
%   WEIGHT  character       Weight type: 'hard' hard-rejection (trimming),
%                                        'soft' soft-rejection (smooth downweighting)
%
%OUTPUT:
%   R       (n x 1)         Radii
%   W       (n x 1)         Weights based on trimmed radii
%   MU      (1 x m x p)     Trimmed mean
%   PC      (q x m x p)     Trimmed principal components
%   L       (q x 1)         Eigenvalues
%   SC      (n x q)         Component scores
%

%% Initialization
if nargin~=6
    error('Incorrect number of input variables')
end
if ~(strcmp(weight,'hard')||strcmp(weight,'soft'))
    error('Variable WEIGHT must be either "hr" or "sr"')
end
[n,m,p] = size(x);
if size(t,2)~=m
    if all(size(t)==[m,1])
        t = t';
    else
        error('Dimensions of T and X inconsistent')
    end
end
h = [(t(2)-t(1))/2,(t(3:m)-t(1:m-2))/2,(t(m)-t(m-1))/2];
na = ceil(alpha*n);
nb = n-floor(beta*n);

%% Interdistances and radii
r = zeros(n,1);
for i = 1:n
    d = sqrt(sum((x - repmat(x(i,:,:),[n,1,1])).^2,3)*h');
    d = sort(d);
    r(i) = d(na);
end

%% Weights
w = zeros(n,1);
d = sort(r);
rnk_n = sum(ones(n,1)*r'<=r*ones(1,n),2)/n;
switch weight
    case 'hard'
        w(r<d(nb)) = 1;
        if beta==0
            w = ones(n,1);
        end
    case 'soft'
        a = 0.50;
        b = 1-beta;
        w(rnk_n<=a) = 1;
        i = find(rnk_n>a&rnk_n<=b);
        w(i) = (rnk_n(i)-b).*((1/(a-b))+(rnk_n(i)-a).*(2*rnk_n(i)-(a+b))/(b-a)^3);
end

%% Estimators
mu = sum(repmat(w,[1,m,p]).*x,1)/sum(w);
pc = [];
l = [];
sc = [];
if q>0
    xw = repmat(sqrt(w),[1,m,p]).*(x-repmat(mu,[n,1,1]));
    G = zeros(n,n);
    for i = 1:n
        G(:,i) = sum(repmat(xw(i,:,:),[n 1 1]).*xw,3)*h'; % Gram matrix
    end
    opts.disp = 0;
    [U,D] = eigs(G,q,'LM',opts);
    l = diag(D);
    pc = zeros(q,m,p);
    sc = zeros(n,q);
    for i = 1:p
        pc(:,:,i) = (U*diag(1./sqrt(l)))'*xw(:,:,i);
        sc = sc + (x(:,:,i)-repmat(mu(:,:,i),[n 1]))*diag(h)*pc(:,:,i)';
    end
    l = l/sum(w);
end