function [alpha,w,b]=reestimate_phi(phi,Y,alpha0, alpha_prune, gamma_a,gamma_b)
% check dimension of inputs
  [n,k]=size(phi{1});
  assert(length(Y)==n);
  
  
% proceed as usual
  alpha=alpha0;
  w=zeros(size(alpha));
  
  numiter=0;
  maxnumiter=10;
  terminate=0;
  eps_tol=1E-5;

  while ~terminate
    tick=cputime;
    numiter=numiter+1;
    alpha_old=alpha;
    inind=(abs(alpha)<alpha_prune);   %% ALPHA may become negative!!
    fprintf('iter=%i\tnum pruned features=%i\n',numiter, sum(~inind));

%    [w,b, invH]=logist_phi(phi,Y,alpha,inind,w);
    [w,b, invH]=logist_phi(phi,Y,alpha,inind);
    dlogdetH=diag(invH);
    
    gamma=1-alpha(inind).*dlogdetH;
    alpha(inind)=(gamma+2*gamma_a)./(w(inind).^2+2*gamma_b);
%    alpha=abs(alpha);  %% CLAMP ALPHA TO BE POSTIIVE!
    alpha(alpha<0)=inf;
    
    if (norm(1./(alpha+eps)-1./(alpha_old+eps))<eps_tol | numiter >= maxnumiter)
      terminate=1;
    end
  end
  

% last run to ensure wb agrees with alpha
  inind=(abs(alpha)<alpha_prune);   %% ALPHA may become negative!!
  [w,b, invH]=logist_phi(phi,Y,alpha,inind,w);

  w(~inind)=0;
