function [r,s_pred,J] = l1_l2_irls(s,w,lambda,Max_iter,tol,silent);
% l1_l2_irls: The l1 regularization term plus l2 misfit cost function is minimized
%             to solve the sparse deconvolution problem.
%
% Find r by minimizing the cost 
%
%        J =   l2-Misfit      + lambda l1-Norm
%          = | w * r - s|_2^2 + lambda |r|_1                  (*)
%
% where w * r is the convolution of a known wavelet (w)  and the unknown reflectivity
% series (r), the observation vector (s) is the seismogram and lambda is the trade-off parameter. 
% 
% IN.
%
%  s(ns):        the seismogram (a column vector)
%  w(nw):        the source wavelet (a column vector) with nw<<ns
%  lambda:       regularization parameter
%  Niter:        maximum number of IRLS iterations 
%  tol:          stop if the relative difference of the l2 norm of the solution between 
%                consecutive iterations is less than tol
%  silent:       1 -> silent
%                0 -> print cost and relative difference of the l2 norm 
%
% OUT.
%
%  r(ns):        the sparse reflectivity
%  s_pred:       the predicted trace, this is the convolution of the
%                wavelet with the estimated reflectivity
%  J:            Matrix where the first dimension is iteration number, second dimension is 
%                J(:,1): cost function (*)
%                J(:,2): l2 misfit
%                J(:,3): l1 regularization term
%                J(:,4): relative diffenrece of the l2 norm of the solution at 
%                        consecutive iterations
%
%  Please refer this article when using the code: 
%
%  M.D. Sacchi, M.D., 1997, Re-weighting strategies in seismic 
%  deconvolution: Geophysical Journal International, 129, 651-656.
%
%  M D Sacchi
%  SAIG  
%  Department of Physics
%  University of Alberta
%
%  (c) M D Sacchi, 2012
%
%
%
%
    if silent == 0; 
     header = '           k           J                 RND';
     format = '       %5.0f  %16.7f  %16.7f';
     disp('  ');      disp(header);
    end


 ns = length(s); 
 nw = length(w); 
 nr = ns;
 ns_padded = ns+nr+nw-1; 
 s = [s;zeros(nw-1,1)];

 W = convmtx(w,nr); 

 g = W'*s;
 R = W'*W;

% Initialization 

  Q = lambda*eye(nr);      k = 1;  
  info = 0;                rnorm_old = 0;
  epsilon = 1.e-5; 

 while (k<=Max_iter) & (info == 0);

   r = (R + lambda*Q)\g;

   rnorm = norm(r,2); 
   delta  = abs(rnorm-rnorm_old)/rnorm;
   rnorm_old = rnorm;
   info = (delta<tol);  

   q = 1./(abs(r)+epsilon);
   Q = diag(q);

   s_pred = conv(w,r);

   Misfit  = sum((s_pred-s).^2) ;
   l1_norm = sum(abs(r)); 

   J(k,1)  = Misfit + lambda*l1_norm;
   J(k,2)  = Misfit;
   J(k,3)  = l1_norm;
   J(k,4)  = delta;

   if silent == 0; 
    disp( sprintf(format, k, J(k,1), J(k,4)));
   end;

   k = k+1;

 end;

 s_pred = s_pred(1:ns); 

 if info==1; disp(sprintf('Relative Norm Difference <= tol                info = %d ',info)); end;
 if info==0; disp(sprintf('k = maximum number of iterations  was reached  info = %d ',info)); end;

return;