% computes the Variational Garrote solution for all gamma
%
%           http://arxiv.org/abs/1109.0486
%
% see train_VG.m
%
% Copyright (C) 2011-, Vicenç Gomez, Bert Kappen
% Last modified: Nov 1, 2013 Vicenç Gomez

min_total_val = inf;
for ini = 1:ninits_vg,    
    %% inits
    minklmin = inf(ng,1);
    v_mf_all = inf(ng,n);
    m_mf_all = inf(ng,n);
    iter_all = inf(ng,2);
    beta_mf_all = inf(1,ng);
    error_mf_all = inf(1,ng);
    errorv_mf_all = inf(1,ng);
    kl_all = inf(ng,2);
    v_all = inf(ng,2,n);
    m_all = inf(ng,2,n);
    beta_all = inf(ng,2);
    imin = inf(1,ng);

    %% random init m
    m = 0.5*ones(1,n) + 0.1*randn(1,n); 
    m(m<=0.2) = 0.2;
    m(m>=0.8) = 0.8;
    eps2 = .4;
    while sum(m)>max_sum_m
        m = eps2*ones(1,n) + 0.1*randn(1,n);
        eps2 = eps2/2;
    end

    %% the estimated inverse noise variance beta is initialized as the output variance
    beta = 1/sigmay;
    
    % change this for more than one fwd-bwd pass
    for numpass = 1,
        
        %% for gamma in gamma_min to gamma_max or when some criteria are satisfied
        i=0;
        while (beta<beta_max) && (i<ng) && (sum(m)<max_sum_m)
            i = i+1; 
            gamma = gamma_all(i)
            tic
            eval(METHOD);
            toc
            v_all(i,1,:) = v;
            m_all(i,1,:) = m;
            beta_all(i,1) = beta;
            iter_all(i,1) = iter;
            kl_all(i,1) = kl1;
            %fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
        end;
        if beta>=beta_max
            fprintf('\tbeta > beta_max (%.3f > %.3f)\n', beta, beta_max);
        end
        if sum(m)>=max_sum_m
            fprintf('\tsum(m) > max_sum_m (%.3f > %.3f)\n', sum(m), max_sum_m);
        end
        if i>1
            m = squeeze(m_all(i-1,1,:))'; 
        end

        %% for gamma is current gamma decreasing to gamma_min 
        imax=i-1;
        for i=imax:-1:1,
            gamma = gamma_all(i);
            eval(METHOD);
            v_all(i,2,:) = v;
            m_all(i,2,:) = m;
            beta_all(i,2) = beta;
            iter_all(i,2) = iter;
            kl_all(i,2) = kl1;
            %fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
        end;
    end
   
    %% take the solution with smallest validation error
    for i=1:imax,
        err_fwd = 1/pv*sum((yv-squeeze(v_all(i,1,:))'*xv).^2,2);
        err_bwd = 1/pv*sum((yv-squeeze(v_all(i,2,:))'*xv).^2,2);
        if err_fwd<err_bwd
            imin(i) = 1;
        else
            imin(i) = 2;
        end
        v_mf_all(i,:) = squeeze(v_all(i,imin(i),:))';
        m_mf_all(i,:) = squeeze(m_all(i,imin(i),:))';
        beta_mf_all(i) = beta_all(i,imin(i));
        error_mf_all(i) = 1/p*sum((y-v_mf_all(i,:)*x).^2,2);
        errorv_mf_all(i) = 1/pv*sum((yv-v_mf_all(i,:)*xv).^2,2);
    end
    [minerrorv i] = min(errorv_mf_all(1:imax));

    %% keep best solution so far in terms of validation error
    fprintf('(%d/%d) ', ini, ninits_vg);
    if minerrorv<min_total_val

        min_total_val = minerrorv;
        gamma_mf = gamma_all(i);
        v_mf = v_mf_all(i,:);
        m_mf = m_mf_all(i,:);
        error_mf = error_mf_all(i);
        errorv_mf = errorv_mf_all(i);
        l1_mf = sum(abs(v_mf));
        if isnan(errorv_mf)
            keyboard;
        end
        beta_mf=beta_mf_all(i);
        if (beta_mf==beta_max)
            fprintf('beta_max too small: beta_mf %6.4f, beta_max %6.4f\n', beta_mf,beta_max);
        end;
        if (gamma_mf==gamma_min)
            fprintf('gamma at minimum range boundary\n');
        end;
        if (gamma_mf==gamma_max)
            fprintf('gamma at maxium range boundary\n');
        end;

        %% plot figures if required
        if plotfigures

            figure;
            subplot(3,2,1)
            plot(gamma_all,kl_all,gamma_mf,0,'*r')
            xlabel('\gamma')
            ylabel('F');
            legend('fwd','bwd','location','southwest')
            subplot(3,2,2)
            plot(gamma_all,sum(m_mf_all,2),gamma_mf,0,'*r')
            xlabel('\gamma')
            ylabel('sum m')
            subplot(3,2,3)
            plot(gamma_all,beta_all,gamma_mf,0,'*r')
            xlabel('\gamma')
            ylabel('\beta')
            subplot(3,2,4)
            plot(gamma_all,error_mf_all,gamma_all,errorv_mf_all,gamma_mf,1e0,'*r')
            ylabel('error')
            legend('train','val','location','northwest')
            xlabel('\gamma')
            set(gca,'yscale','log');
            drawnow

            if toy,
                [ws, idx] = sort(abs(w),'ascend');
            else
                idx = 1:n;
            end
            
            subplot(3,2,5);
            first_nzero = find(m(idx)>.1,1,'first');
            if isempty(first_nzero)
                first_nzero = 1;
            end
            if toy
                bar([m_mf(idx)' (w(idx)~=0)']);
                legend('m', 's','location','northwest');
            else
                bar(m_mf(idx));
                legend('m','location','northwest');
            end
            set(gca,'xticklabel',[]);
            set(gca,'ylim',[0 1]);
            hold on;
            plot([0 n],[.5 .5],'--k');
            set(gca,'xlim',[first_nzero n]);
            subplot(3,2,6);
            if toy
                bar([v_mf(idx)' w(idx)']);
                legend('v', 'w^*','location','northwest');
            else
                bar(v_mf(idx));
                legend('v','location','northwest');
            end
            set(gca,'xlim',[first_nzero n]);
            set(gca,'xticklabel',[]);
            drawnow;
 
            if toy
                %% compute ROCs
                w2 = w>0;
                [tpr,fpr,thresholds] = roc(w2,abs(v_mf));
                tp = [0 ; tpr(:) ; 1];
                fp = [0 ; fpr(:) ; 1];
                nn = size(tp, 1);
                auc_v_mf = sum((fp(2:nn) - fp(1:nn-1)).*(tp(2:nn)+tp(1:nn-1)))/2;
                fprintf('auc = %.3f\t', auc_v_mf);
                
                fprintf('L1 err = %.3f\t', sum(abs(v_mf-w)));
                fprintf('beta = %.3f', beta_mf_all(i));
            end
        end
    end
    fprintf('\n');
    if imax<0
        m_mf = m;
        v_mf = v;
        gamma_mf = gamma;
    end
end
