function [ outW, outNumIter, outObj, out ] = L2RtrR21_W( inX, inY, inPara, inW0)
% reweighted method to solve the following problem
% ||Y- X'W||_F^2 + lambda1*||W||_* + lambda2*||W||_21
% input: 
%       inX: data :d x n  
%       inY: label:n x m
%       inPara: 
%               inPara.lambda1: regularization  parameter for the trace
%               norm
%               inPara.lambda2: regularization  parameter for the 2,1 norm
%               inPara.maxIter: the max number of iteration
%               inPara.thresh: stopping criterion threshold
%       inW0: initialized weight matrix d x c (optional)
% output:
%       outW: the weigth matrix : d x m
%       outObj: the output obj
%       outNumIter: the number of iteration untill converge
%       out: structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check the size of the input arguments
%inX = inX(1:100, :);
[d, n1] = size(inX);
[n2, m] = size(inY);
if(n1 ~= n2)
    error('The size of the input arguments are not compatible !');
else
    n = n1;
    clear n1 n2
end
% parameter settings
maxIter = inPara.maxIter;
thresh = inPara.thresh;
lambda1 = inPara.lambda1;
lambda2 = inPara.lambda2;
% initialization
if nargin < 4
    D1 = diag(ones(d,1));
    D2 = diag(ones(d,1));
else
    W = inW0; % if there is W initialization in the input argument
    D1 = 0.5*(W*W')^(-0.5);
    Wi2 = sqrt(sum(W.*W, 2) + eps);
    d2 = 0.5./(Wi2);
    D2 = diag(d2);
end

obj = zeros(maxIter, 1);
out.loss = zeros(maxIter, 1);
out.termTr = zeros(maxIter, 1);
out.term21 = zeros(maxIter, 1);
% loop Alg
for t = 1: maxIter
    fprintf('processing iteration %d...\n', t);
    % fix D1 and D2, update W
    W = (inX*inX' + lambda1*D1 + lambda2*D2)\(inX*inY);
    W = real(W);
    % fix W, update D1
    %D1 = 0.5*(W*W')^(-0.5);  %% may cause D1 not real number %%
	[U S V] = svd(W); s = diag(S);%s = zeros(n,1); t = diag(S); s(1:length(t)) = t;
    D1 = 0.5*U*diag((s.*s).^(-0.5))*U';
    % fix W, update D2
    Wi2 = sqrt(sum(W.*W, 2) + eps);
    d2 = 0.5./(Wi2);
    D2 = diag(d2);
    D2 = sparse(D2);
    % calculate trace Norm
    temp = W'*W;
    temp1 = temp.^0.5;
    TrW = trace(temp1);
    % calculate 21-Norm
    W21 = sum(Wi2);
    % calculate obj
    temp2 = norm(inY - inX'*W, 'fro');
    out.loss(t) = temp2^2;
    out.termTr(t) = lambda1 * TrW;
    out.term21(t) = lambda2 * W21;
    obj(t) = out.loss(t) + out.termTr(t) + out.term21(t);
    % stopping criterion
    if(t > 1)
        diff = obj(t-1) - obj(t);
        if(diff < thresh)
            break;
        end
    end
end
outW = W;
% outW = outW.*(abs(outW) > thresh);
outNumIter = t;
outObj = obj(1:outNumIter);
out.loss = out.loss(1:outNumIter);
out.termTr = out.termTr(1:outNumIter);
out.term21 = out.term21(1:outNumIter);


% debug
%figure(100), plot(obj(1:outNumIter));
end


