function [L,S,Out] = TRPCA_WSTNN_ADMM(X,opts)

% Solve the WSTNN-based TRPCA problem by ADMM

if ~exist('opts', 'var')
    opts = [];
end    
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'max_rho');      max_rho = opts.max_rho;        end
if isfield(opts, 'beta');        beta = opts.beta;                end
if isfield(opts, 'max_beta');      max_beta = opts.max_beta;        end
if isfield(opts, 'alpha');        alpha = opts.alpha;                end
if isfield(opts, 'gamma');         gamma = opts.gamma;              end
if isfield(opts, 'lambda');        lambda = opts.lambda;                end
Nway = size(X);
%% Initialization

N = ndims(X);
L = zeros(Nway);

Z = cell(N,N);
for i=1:N
    for j=1:N
    Z{i,j} = zeros(Nway); %% the auxiliary tensor
    end
end
P = Z;
S = zeros(Nway);
M = zeros(Nway);

temp = Z;

Out.Res=[]; Out.PSNR=[];
for iter = 1 : max_iter
    %% Let

    Lold = L;
    
    %% update Z
    tau = alpha./beta;
    for i=1:N-1
        for j=i+1:N
            tempL   = tensor_permute(L,Nway,i,j);
            tempP   = tensor_permute(P{i,j},Nway,i,j);
            Z{i,j}  = tensor_ipermute(prox_tnn_my(tempL + tempP/beta(i,j),tau(i,j)),Nway,i,j);
            temp{i,j} = Z{i,j}-P{i,j}/beta(i,j);
        end
    end
    
    
    %% update L
    tempsum = zeros(Nway);
    for i=1:N-1
        for j=1+i:N
        tempsum = tempsum+ beta(i,j)*temp{i,j};
        end
    end
    L = (tempsum+rho*(X-S+M/rho))/(rho+sum(beta(:)));
        
    %% update S
    S = prox_l1(X-L+M/rho,lambda/rho);
  
    %% check the convergence
    dM = X-L-S;
    chg =norm(L(:)-Lold(:))/norm(Lold(:));
    
    Out.Res = [Out.Res,chg];
    if isfield(opts, 'Xtrue')
        XT=opts.Xtrue;
        psnr=PSNR(L,XT);
        Out.PSNR = [Out.PSNR,psnr];
    end
    
    if iter == 1 || mod(iter, 10) == 0
         if isfield(opts, 'Xtrue')
            fprintf('WSTNN: iter = %d   PSNR=%f   res=%f  \n', iter, psnr, chg);
         else
            fprintf('WSTNN: iter = %d   res=%f  \n', iter, chg);
         end       
    end
    
    if chg < tol
        break;
    end
    %% update M & P
    for i=1:N-1
        for j=i:N 
           P{i,j} = P{i,j}+beta(i,j) * (L-Z{i,j});
        end
    end
    M = M + rho*dM;
    beta = min(gamma*beta,max_beta);    
    rho = min(gamma*rho,max_rho);  
end