function [T,C,w] = GMt(U,V,lambdas,nu,weight,alpha,dsp)
%function [T,C,w] = GMt(U,V,l,nu,wgt,trp,dsp)
% GMt functional regression estimator
%
% INPUT
%   U   (n x p)     Reduced-rank standardized component scores of X
%   V   (n x p)     Reduced-rank standardized component scores of Y
%   l   (p x 1)     Eigenvalues of PCs of X
%   nu  (scalar)    Degrees of freedom of the t used as loss function
%   wgt (char)      Trimming type: "metric" or "rank"
%   trp (scalar)    Trimming proportion (between 0 and 1)
%   dsp (char)      Display of iteration progress: "on" or "off"
%
% OUTPUT
%   T   (p x q)     Regression estimator
%   C   (pq x pq)   Asymptotic covariance matrix of sqrt(n)*vec(T)
%   w   (n x 1)     Individual weights

% Initialization
[n,p] = size(U);
q = size(V,2);
rho = @(x)((nu+q)*log(1+x/nu));
psi = @(x)((1+q/nu)./(1+x/nu));

% Weight computation
lmb1 = zeros(p,1);
i = find(lambdas>0);
lmb1(i) = 1./sqrt(lambdas(i));
maha = sum((U*diag(lmb1)).^2,2);
switch lower(weight)
    case 'metric'
        w = maha<=chi2inv(1-alpha,p);
    case 'rank'
        rnk = sum((ones(n,1)*maha')<=(maha*ones(1,n)),2);
        w = (rnk/n)<=(1-alpha);
    otherwise
        error('WEIGHT must be "metric" or "rank"')
end

% Iterative computation of estimators
T = zeros(p,q);
R = V-U*T;
S = diag(median(R.^2));
D = w.*sum((R/S).*R,2);
F = n*log(det(S)) + sum(rho(D));
if strcmp(dsp,'on')
    disp(['Initial obj. func. F = ' num2str(F)])
end
err = 1;
iter = 0;
while err>1e-3 && iter<100
    iter = iter + 1;    
    F0 = F;
    v = psi(D).*w;
    A = U'*(repmat(v,[1 p]).*U);
    B = U'*(repmat(v,[1 q]).*V);
    T = A\B;
    R = V-U*T;
    S = R'*(repmat(v,[1 q]).*R)/n;
    D = w.*sum((R/S).*R,2);
    F = n*log(det(S)) + sum(rho(D));
    err = abs(F-F0)/F0;
    if strcmp(dsp,'on')
        disp(['Iter ' num2str(iter) ', F = ' num2str(F) ', Error: ' num2str(err)])
    end
end

% Asymptotic covariance of regression estimator ("sandwich" formula)
if nargout>=2
    psi1 = @(x)(-(1/nu+q/nu^2)./(1+x/nu).^2);
    v = psi(D).*w;
    v1 = psi1(D).*w.^2;
    A = U'*(repmat(v,[1 p]).*U)/n;
    BB = zeros(p*q);
    CC = zeros(p*q);
    for i = 1:n
        BB = BB + v(i)^2*kron(R(i,:)'*R(i,:),U(i,:)'*U(i,:));
        CC = CC + v1(i)*kron(R(i,:)'*R(i,:),U(i,:)'*U(i,:));
    end
    BB = BB/n;
    CC = CC/n;
    AA = 2*kron(inv(S),eye(p))*CC + kron(eye(q),A);
    C = (AA\BB)/AA;
end