function [labels ,converge_Z ,converge_Z_G,Z,Zbar,U]=STONE(X,gt,A,alpha,beta,gamma,delta)
numC = length(unique(gt));
V = length(X); %number of views c=size(A,2) 璇ヨ鍙ヨ繑鍥炵殑鏄煩闃礎鐨勫垪鏁?
N = size(X{1},2);% number of data points
[~,M]=size(A{1});% dv by M
mu1 = 10e-5; max_mu1 = 10e10; 
mu2 = 10e-5; max_mu2 = 10e10; 
rho = 10e-5; max_rho = 10e10; 
pho = 2;
% rho = 0.1;
sX = [N,M,V];
%Initilize Z,Ls,tensor J,multiplier W
for i = 1:V
    B{i}=zeros(N,M);
    Lh{i}=zeros(M,M);
    dv=size(X{i},1);
    Z{i} = zeros(N,M);
    G{i} = zeros(dv,dv);
    J{i} = zeros(N,M);
    K{i} = zeros(N,M);  %multiplier
    Q{i} = zeros(N,M);  %multiplier
    R{i} = zeros(dv,N);
    E{i} = zeros(dv,N);
end

% NITER = 30;
%
% for iter=1:NITER
iter = 0;Isconverg = 0; epson = 1e-7;
start=1;
converge_Z=[];converge_G=[];converge_Z_G=[];
while(Isconverg == 0)
    fprintf('----processing iter %d--------\n', iter+1);
    %% ============================== Upadate Lh{i} ===============================
    for i = 1 :V
        if start==1
          H{i} = constructW_PKN(A{i}, 3);
          Dh = diag(sum(H{i},1)+eps);
          Lh{i} = eye(M,M)-Dh^-0.5*H{i}*Dh^-0.5;
        else
        %------------modified to hyper-graph---------------
          param.k = 3;
          HG = gsp_nn_hypergraph(A{i}', param);
          Lh{i} = HG.L;
        end
        
    end
    start=0;
    %% =========================== Upadate Z{i} ===========================
    Zbar=[];
    for i=1:V
        Z{i} = (rho*J{i}-K{i}+mu1*X{i}'*A{i}-mu1*X{i}'*G{i}'*A{i}-mu1*E{i}'*A{i}+R{i}'*A{i}+mu2*Q{i}-B{i})/((rho+mu2)*eye(M)+mu1*A{i}'*A{i});
        Zbar=cat(2,Zbar,1/sqrt(V)*Z{i});
    end
    
    %% =========================== Upadate G{i} ===========================
    for i=1:V
        dv=size(X{i},1);
        temp=X{i}-A{i}*Z{i}'-E{i}+R{i}/mu1;
%         G{i}=(mu1*temp*X{i}')/(alpha*eye(dv)+mu1*X{i}*X{i}');
        G{i}=(mu1*temp*X{i}')*pinv((alpha*eye(dv)+mu1*X{i}*X{i}'));
    end
    
    %% ====================== Upadate E{i} ======================
    D = [];
    for i=1:V
        D = [D;X{i}-A{i}*Z{i}'-G{i}*X{i}+R{i}/mu1];
    end
    [Econcat] = solve_l1l2(D,beta/mu1);
    
    beg_ind = 0;
    end_ind = 0;
    for i=1:V
        if(i>1)
            beg_ind = beg_ind+size(X{i-1},1);
        else
            beg_ind = 1;
        end
        end_ind = end_ind+size(X{i},1);
        E{i} = Econcat(beg_ind:end_ind,:);
    end
    %% ====================== Upadate Q{i} ======================
    for i=1:V
        Q{i}=(mu2*Z{i}-B{i})/(2*gamma*Lh{i}+mu2*eye(M));
    end

    %% =========================== Upadate A{i} ===========================
    for i=1:V
        tempB{i}=mu1/2*(X{i}-G{i}*X{i}-E{i}+R{i}/mu1)*Z{i};
        [U,~,Vb]=svd(tempB{i},'econ');
        A{i}=U*Vb';
    end
    %% =========================== Upadate J ===========================
    Z_tensor = cat(3, Z{:,:});%C = cat(dim, Z, B) �?? dim 鎸囧畾鐨勬暟缁勭淮搴︿覆鑱旀暟缁? Z �?? B
    K_tensor = cat(3, K{:,:});
    z = Z_tensor(:);
    k = K_tensor(:);
    %     [j, ~] = wshrinkObj(a+1/rho*o,beta/rho,sX,0,3,omega);
    %     [j, ~] = wshrinkObj_weight_lp(a+1/rho*o,beta/rho*omega,sX,0,3,p);
    %     J_tensor = reshape(j, sX);
    %      [J_tensor,objV]=solve_J(a+1/rho*o,rho/beta,sX,delta);
    [J_tensor,objV]=solve_J(z+1/rho*k,rho,sX,delta);
    for i=1:V
        J{i} = J_tensor(:,:,i);
    end

    
    %% =========================== Upadate multiplier ===========================
    for i=1:V
        K{i}=K{i}+rho*(Z{i}-J{i});
        R{i}=R{i}+mu1*(X{i}-A{i}*Z{i}'-G{i}*X{i}-E{i});
        B{i}=B{i}+mu2*(Z{i}-Q{i});
    end
    mu1 = min(mu1*pho, max_mu1);
    mu2 = min(mu2*pho, max_mu2);
    rho = min(rho*pho, max_rho);
    %% ===========================Checking Coverge Condition ===========================
    
    Isconverg = 1;
    sum_converge_Z=0;
    sum_converge_G=0;
    sum_converge_Z_G=0;
    for i=1:V
        if(norm(Z{i}-J{i},inf)>epson)
            history.norm_Z=norm(Z{i}-J{i},inf);
            fprintf('    norm_Z %7.10f   \n ', history.norm_Z);
            Isconverg = 0;
            sum_converge_Z=sum_converge_Z+history.norm_Z;
        end
%         if(norm(G{i}-K{i},inf)>epson)
%             history.norm_G=norm(G{i}-K{i},inf);
%             fprintf('    norm_G %7.10f   \n ', history.norm_G);
%             Isconverg = 0;
%             sum_converge_G=sum_converge_G+history.norm_G;
%         end
        if(norm(X{i}-A{i}*Z{i}'-G{i}*X{i}-E{i},inf)>epson)
            history.norm_Z_G=norm(X{i}-A{i}*Z{i}'-G{i}*X{i}-E{i},inf);
            fprintf('    norm_Z_G %7.10f   \n ', history.norm_Z_G);
            Isconverg = 0;
            sum_converge_Z_G=sum_converge_Z_G+history.norm_Z_G;
        end
    end
    
    
    converge_Z=[converge_Z sum_converge_Z];
%     converge_G=[converge_G,sum_converge_G];
    converge_Z_G=[converge_Z_G,sum_converge_Z_G];
    
    if sum_converge_Z==inf  ||sum_converge_Z_G==inf
        break
    end
    if (iter>50)
        Isconverg  = 1;
    end
    iter = iter + 1;
    
    
end
%% =========================== get result ===========================
[U,~,~] = mySVD(Zbar,numC);
rng(5489,'twister');
labels=litekmeans(U, numC, 'MaxIter', 100,'Replicates',10);