% dual method for the Variational Garrote
%
% Copyright (C) 2011-, Vicenç Gomez, Bert Kappen
% Last modified: Nov 1, 2013 Vicenç Gomez

iter = 1;
dm = 1;
eta = eta0;
eta_w = eta_w0;
m = m+1e-13*randn(1,n);

precism =1e-13;
m=max(m,precism);
m=min(m,1-precism);
final = false;

kl_tmp=inf(1,maxiter);
if (beta>=beta_max), beta=beta_max-1; end
while ((dm>dmmin)&&(iter<maxiter)&&(beta<beta_max)&&~final) 
	A=eye(p)+1/p*x'*spdiags(m'./(1-m)'./chi_ii,0,n,n)*x;
    klad=A\y';
	beta=p/(y*klad);
	lambda=beta*klad;
   	mold=m;
	klad1=x*lambda;
   	m=(1-eta)*m+eta*sigmoid1(gamma+(klad1'.^2)./chi_ii'./((1-m).^2)/2/beta/p);
	m=max(m,precism);
	m=min(m,1-precism);
   	dm=max(abs(m-mold));

    kl_tmp(iter)=-p/2*log(beta)-1/2/beta*lambda'*A*lambda+y*lambda - gamma*sum(m)+sum(m.*log(m)+(1-m).*log(1-m));
    if iter>1 
        dkl =(kl_tmp(iter)-kl_tmp(iter-1));
        while dkl>precism
            eta=eta/2;
            if eta < precism
                %fprintf('\t[%d] dkl=%e eta=%.3e, dm=%.3e\n',iter,kl_tmp(iter)-kl_tmp(iter-1),eta,dm);
                final=true;
                m=mold;
                iter=iter-1;
                break;
            end
            m=(1-eta)*m+eta*sigmoid1(gamma+(klad1'.^2)./chi_ii'./((1-m).^2)/2/beta/p);
            m=max(m,precism);
            m=min(m,1-precism);
            dm=max(abs(m-mold));
            
            kl_tmp(iter)=-p/2*log(beta)-1/2/beta*lambda'*A*lambda+y*lambda - gamma*sum(m)+sum(m.*log(m)+(1-m).*log(1-m));
        end
    end
    if ~final
        iter=iter+1;
    end
end
if iter==maxiter%, gamma, iter, 
    disp('no convergence');
end
kl1=-p/2*log(beta)-1/2/beta*lambda'*A*lambda+y*lambda - gamma*sum(m)+sum(m.*log(m)+(1-m).*log(1-m));
v=1./chi_ii'.*m./(1-m).*klad1'/beta/p;
