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

%   DANSER tries to solve the following problem

%   sum (1/2)*(sum_l=1^L||y_l-Hs_l||_2)^p + lambda sum ||s^i||_2^p + gamma ||H -Z||_F^2
%   st   ||z_i-b_i||<=delta

% where the inputs are :
% B: a given dictionary, possibly with erros
% H: initialization of a sought "clean" dictionary (initialization H=B)
% S: initialization of the abundance map
% p: 0<p<1; recommended: p=0.5, 0.75.
% lambda: a large lamabda promotes row-sparser solution of S
% gamma: typically gamma = 1000 or 10000.
% Y: an M x L hyperspectral image, where M is the number of spectral bands,
% and L is the number of pixels.

% reference:
%  X. Fu, W.-K. Ma, J. M. Bioucas-Dias, and T.-H. Chan,
% ''Semiblind Hyperspectral Unmixing in the Presence of Spectral Library Mismatches,''
%  IEEE Transactions on Geoscience and Remote Sensing, to appear, 2016.
%  coded by Xiao Fu email: xfu@umn.edu


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


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;zeros(K,L)]'*H_tilde;
    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; % here the stopping criterion can be changed to accomandate your problem
        break
    end
    
    S_pre = S;
    
end


end