function [r,w,mu,pc,l,sc] = fte2(s,t,x,q,alpha,beta,weight)
%function [r,w,mu,pc,l,sc] = fte2(s,t,x,q,alpha,beta,weight)
%Functional Trimmed Estimators (for bivariate, real-valued functions)
%
%INPUT:
%   S       (1 x m1)        First input variable
%   T       (1 x m2)        Second input variable
%   X       (m1 x m2 x n)   Data surfaces
%                           (To accomodate non-rectangular domains, NaN's are accepted
%                            but are replaced by zero - so do NOT use NaN's for missing values)
%   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      (m1 x m2)       Trimmed mean
%   PC      (m1 x m2 x q)   Trimmed principal components
%   L       (q x 1)         Eigenvalues
%   SC      (n x q)         Component scores
%

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

%% Interdistances and radii
r = zeros(n,1);
d = zeros(n,1);
I = isnan(x);
x0 = x;
x0(I) = 0;
for i = 1:n
    for j = 1:n
        d(j) = sqrt(hs*((x0(:,:,j)-x0(:,:,i)).^2)*ht');
    end
    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
w = reshape(w,[1 1 n]);
mu = sum(repmat(w,[m1 m2 1]).*x,3)/sum(w,3);
mu0 = mu;
I = isnan(mu);
mu0(I) = 0;

pc = [];
l = [];
sc = [];
if q>0
    xw = repmat(sqrt(w),[m1 m1 1]).*(x0-repmat(mu0,[1 1 n]));
    G = zeros(n,n);
    for i = 1:n
        for j = 1:n
            G(i,j) = hs*(xw(:,:,i).*xw(:,:,j))*ht'; % Gram matrix
        end
    end
    opts.disp = 0;
    [U,D] = eigs(G,q,'LM',opts);
    l = diag(D);
    pc = zeros(m1,m2,q);
    sc = zeros(n,q);
    for k = 1:q
        pc_aux = zeros(m1,m2);
        for i = 1:n
            pc_aux = pc_aux + U(i,k)*xw(:,:,i);
        end
        pc_aux = pc_aux/sqrt(l(k));
        for i = 1:n
            sc(i,k) = hs*((x0(:,:,i)-mu0).*pc_aux)*ht';
        end
        pc_aux(I) = NaN;
        pc(:,:,k) = pc_aux;
    end
    l = l/sum(w,3);
end