function [ S, it] = DNANSER_mod(B,H,S,p,Y,lambda,gamma,delta,maxiter)

% Input
% B: spectral dictionary;
% H: the denoised dictionary;
% S: coefficient matrix; row sparse
% Y: data
% p: 0<p<=1, for sparsity promoting, (e.g. p=0.5)
% lambda: regularization parameter, (e.g. lambda = 0.1)
% gamma: reularization paramter (a large one) (e.g., gamma=1000);
% delta: the paramter indicating how bad the dictionary is; larger delta
% means we have more freedom to adjust the dictionary (e.g.,delat=1e-1)

% output
% S: estimated coefficient matrix

%   DANSER tries to solve the following problem
%   sum (1/2)*(||Y-HS||_F^2 + lambda sum ||s^i||_2^p + gamma ||H -Z||_F^2
%   st   ||z_i-b_i||<=delta

[M,K]=size(H);  % This K might be over-estimated
[M,L]=size(Y);

%% initialization
Z = H;
[K,L]=size(S);
eta = 1e-5;
ST = S';
Wt1 =(sum(ST.^2));
Wt = (p/2)*(Wt1+eta).^((p-2)/2);
S_pre = S;

for it=1:maxiter
 
    %------------------update S ---------------------

    H_tilde=[H;sqrt(2*lambda.*diag(Wt))];
     
    W = Y'*H;
    V =  H_tilde'*H_tilde;
    
    for k=1:K
        s= S(k,:).'+(W(:,k) - (S.')*V(:,k))/V(k,k);
        S(k,:)=max(s.',0);
    end

    %------------------update W (closed form)---------------------
    ST = S';
    Wt1 =(sum(ST.^2));
    Wt = (p/2)*(Wt1+eta).^((p-2)/2);
 
    %------------------update H (closed form) ----------------------
    H = (Y*S' + gamma*Z)/(S*S'+ gamma*eye(K));

    
    % --------------   update Z (closed-form)-------------------
    z_h = sqrt(sum((H-B).^2));
    ind_1=find(z_h<=delta);
    Z(:,ind_1)=H(:,ind_1);
    ind_2=find(z_h>delta);
    Z(:,ind_2) = B(:,ind_2)  +  delta*(H(:,ind_2)- B(:,ind_2))*diag(1./z_h(ind_2));
    
    
    
    %--------------- update the objective function --------------------
    relative_err = sum(sum((S - S_pre).^2));
    
    if relative_err < 1e-5
        break
    end
    S_pre = S;
end

end
