
% lsdt.m function for computing sparse coefficients and subspace mapping.

function [C,A,obj,K,Ks,Kt,Kt_test]=lsdt(Xt,Xs,X,Xt_test,dim,lambda1,lambda2,maxiter,Kernel,kervar1,kervar2)


numtrain    = size(X,2) ;
K           = gram(X',X',Kernel,kervar1,kervar2);
Kt           = gram(X',Xt',Kernel,kervar1,kervar2);
Ks           = gram(X',Xs',Kernel,kervar1,kervar2);
Kt_test      = gram(X',Xt_test',Kernel,kervar1,kervar2);
% K           = normalizeK(K);
% Ks=K(:,1:size(Xs,2));
% Kt=K(:,size(Xs,2)+1:size(Xs,2)+size(Xt,2));
% Kt_test=K(:,size(Xs,2)+size(Xt,2)+1:end);
[A Beta]    = eig(K) ;
[dummy id]  = sort(diag(Beta),'descend');
A           = A(:,id(1:dim)) ;A=real(A);

[V S dummy] = svd(K) ;
sqrtS       = sqrt(S) ;
inv_sqrtS   = pinv(sqrtS) ;
I           = eye(numtrain);
affine      = true;%false;
alpha       = lambda1;%1e0;
thrr        = 10^(-6);

% start learning
for loop=1:maxiter
    
    %solving for C
%     Z           = A'*K ;
    Zt= A'*Kt;% Kt
    Zs=A'*Ks; % Ks
    
    C       = admmLasso_mat_func(Zt,Zs,affine,alpha,thrr); % LSDT: sparse constraint
     
   % solving for A
    e=1;
%     switch Kernel
%         case 'linear'
% %            M  = lambda1*(X'/(X*X'+e)*Xt-X'/(X*X'+e)*Xs*C)*(X'/(X*X'+e)*Xt-X'/(X*X'+e)*Xs*C)' - lambda2*I ;% linear
%            M  = lambda1*(K\Kt-K\Ks*C)*(K\Kt-K\Ks*C)' - lambda2*I ;
%         case 'gauss'
%            M  = lambda1*(K\Kt-K\Ks*C)*(K\Kt-K\Ks*C)' - lambda2*I ;
% %            M  = lambda1*(X'/(X*X'+e)*Xt-X'/(X*X'+e)*Xs*C)*(X'/(X*X'+e)*Xt-X'/(X*X'+e)*Xs*C)' - lambda2*I ;% linear
%     end % nonlin

%     tmp=inv(K);
%     if tmp(1,1)>10^6
%        M  = lambda1*(pinv(K)*Kt-pinv(K)*Ks*C)*(pinv(K)*Kt-pinv(K)*Ks*C)' - lambda2*I ;
%     else 
%        M  = lambda1*(K\Kt-K\Ks*C)*(K\Kt-K\Ks*C)' - lambda2*I ;
%     end

    M  = lambda1*(pinv(K)*Kt-pinv(K)*Ks*C)*(pinv(K)*Kt-pinv(K)*Ks*C)'-lambda2*I;
    H           = sqrtS*V'*M*V*sqrtS;
    H           = (H+H')/2;
    [G Beta]    = eig(H);%eig(H,H,'chol'); 
    
    %{
    if (any(imag(G(:))) ) %%%% check error
        fprintf(['Complex Eigenvalues' '\n']);
        pause;
    end
    %}
    
    temp        = diag(Beta) ;
    [dummy id]  = sort(temp,'ascend') ;
    G = G(:,id(1:dim)) ; % select eigvectors of smallest eigvalues
    A           = V*inv_sqrtS*G ;
    
    switch Kernel
        case 'linear'
           P=A'*X';obj(loop)=max(diag(C'*C))+lambda1*norm(P*Xt-P*Xs*C,'fro')+lambda2*norm(X-P'*P*X,'fro');
        case 'gauss'
            P=A'*K';obj(loop)=max(diag(C'*C))+lambda1*norm(P*Kt-P*Ks*C,'fro')+lambda2*norm(K-P'*P*K,'fro');
    end
end
