function model = train_VG(X, Y, valp)

%% model = train_VG(X, Y, valp)
%
%           Sparse linear regression using the Variational Garrote method
%           http://arxiv.org/abs/1109.0486
%           
% use train_VG([]) for running with toy data (see make_toy_data.m)
%
% X       : n x p matrix (n number of features, p number of samples)
% Y       : p x 1 vector (outputs)
% valp    : proportion of dataset used for validation (default 0.1)
%
% Copyright (C) 2011-, Vicenç Gomez (v.gomez@science.ru.nl), Bert Kappen
% Last modified: Nov 1, 2013 Vicenç Gomez

    %% parameters
    
    plotfigures = true;     % plot results?
    METHOD = 'dual';%  'regression';%      % method to execute 'dual' or 'regression'

    if isempty(X)
        %% generate example data
        fprintf('Generating toy data ...');
        make_toy_data;
        ok = ~any(isnan(x(:))) && ~any(isnan(xv(:)));
        toy = true;
        fprintf('done\n');
    else        
        %% split train - validation
        if nargin == 2
            valp = .1;
        end
        n = size(X,1);
        p = floor((1-valp)*size(X,2));
        pv = size(X,2)-p;
        idxp = randperm(size(X,2));
        x = X(:,idxp(1:p));
        xv = X(:,idxp(p+1:end));
        y = Y(idxp(1:p));
        yv = Y(idxp(p+1:end));

        %% normalize training data
        x = x-mean(x,2)*ones(1,p);
        dx = sqrt(1/p*sum(x.^2,2));
        x = x./(dx*ones(1,p));
        y = y-mean(y);

        %% normalize validation data
        xv = xv-mean(xv,2)*ones(1,pv);
        dxv = sqrt(1/pv*sum(xv.^2,2));
        xv = xv./(dxv*ones(1,pv));
        yv = yv-mean(yv);

        ok = ~any(isnan(x(:))) && ~any(isnan(xv(:)));
        toy = false;
    end
    sigmay = y*y'/p;  % sigma is output variance
    
    %% other parameters (general)

    epsilon = 1e-4;         % determines gamma_min (smallest m, see Eq. 13)
    ng = 3;                % number of gammas (equivalent to step-size)
    eps_gamma = 1e-4;       % determines gamma_max
    ninits_vg = 1;          % number of random initializations
    dmmin = 1e-10;          % convergence of mean field error 1e-12
    eta0 = 1;               % initial step size of mean field update
    eta_w0 = 0.02;          % initial step size for change in w in dual.m
    maxiter = 1e3;          % maximal number of iterations 1e4
    max_sum_m = n;          % increases gamma until sum(m) == max_sum_m
    beta_max = sigmay;      % increases gamma until beta == beta_max

    if ok
        disp(['method is ' METHOD]);
        if strcmp(METHOD, 'regression')
            if n>2e3	
                fprintf('Computing input-input covariance ...');
            end
            C = x*x'/p;   % C is input data covariance matrix.
            if n>2e3	
                fprintf('done\n');
            end
        end
        b = x*y'/p;       % b is input output covariance
        if pv>0,
            bv = xv*yv'/pv;
            sigmayv = yv*yv'/pv;
        end;
        chi_ii=1/p*sum(x.^2,2); % input data variance
        if sum(abs(chi_ii-1)>1e-10),
            fprintf('input design matrix is not normalized\n');
        end;

        %% set gamma range (min, max and step size)
        gamma_min = log(epsilon*sigmay/p/max(abs(b))); %gamma_min = -(p*max(abs(b))/sigmay) + log(epsilon/(1-epsilon));%
        gamma_max = eps_gamma*gamma_min;
        gamma_all =linspace(gamma_min, gamma_max, ng);


        %% compute VG solution for range of gamma as given in gamma_all
        % first from gamma_min to gamma_max and then in a second pass
        % from gamma_max to gamma_min.
        tic;	
        vg_learn;
        toc

        model.m = m_mf;
        model.v = v_mf;
        model.beta = beta_mf;
        model.gamma = gamma_mf;
        model.x = x;
        model.xv = xv;
        model.y = y;
        model.yv = yv;
        if toy
            model.auc_v_mf = auc_v_mf;
            model.l1err = sum(abs(v_mf-w));
        end
    else
        model = [];
    end

end
