function [retTP,retM] = ValidateTP(trY,trX,vdY,vdX,inFamily,inLink,inOption)
%function [retTP,retM] = ValidateTP(trY,trX,vdY,vdX,inFamily,inLink,inOption)
%
% Determine the tuning parameter using pre-specified validation data.
%
% Model: g(E[Y]) = mu + X*beta;
%
% Input:
%   trY: Y in training data
%   trX: X in training data
%   vdY: Y in validation data
%   vdX: x in validation data
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   inOption:
%       tol -- tolerance,
%       maxit -- maximum number of iterations,
%       bctype -- type of bias correction: 0 ~ no bias correction,
%                 1 ~ approximate (default), type 2 -- optimal;
%       maxcmp -- maximum number of components.
%       maxvar -- maximum number of selected variables;
%       W -- weight matrix;
%       ptype -- type of penalization,
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation,
%                'l1' ~ L1,
%                'scad' ~ SCAD (with gamma=3.7),
%                'mcp' ~ MCP (with gamma=3.7).
%       gamma -- parameter for SCAD or MCP (=3.7 by default);
%       delta -- pace for choosing value of the tuning parameter
%       vc -- cross-validation criterion: 'press' (default),'mr','loglik';
%
% Output:
%   retTP: optimal value of the tuning parameter;
%   retM: The fitting result of the optimal model.
%
% Dabao Zhang     July 1, 2013
%
% Revised by:
%       Dabao Zhang, Oct. 29, 2013 -- Change the initialization of W.
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(5, 7, nargin));

if( nargin<7 ) inOption = []; end
if( ~isfield(inOption,'delta') ) inOption.delta = 0.1; end
if( ~isfield(inOption,'lambda0') ) inOption.lambda0 = 1; end
if( ~isfield(inOption,'vc') ) inOption.vc = 'press'; end
if( ~isfield(inOption,'W') ) inOption.W = []; end
if( (nargin<6)||isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'multibc'
            inLink = 'bclogit';
        case 'multipo'
            inLink = 'pologit';
        case 'multicr'
            inLink = 'crlogit';
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        otherwise
            error('  OptimalTP: The family is not supported...');
    end
end

% Special for multinomial models
bMultiN = 0;
if( sum(strcmp(inFamily,{'multibc','multipo','multicr'})) )
    bMultiN = 1;
end

% Explore all tuning parameter values such that 0<#[beta~=0]<maxvar
idxM = 0;           lambda = inOption.lambda0;
if( isempty(inOption.W) )
    bInitW = 1;
    inOption.W = ones(size(trX,1),1)/size(trX,1);
else
    bInitW = 0;
end
while( idxM<=inOption.maxit )
    disp(['  Tuning Parameter: ' num2str(lambda(end))]);
    idxM = idxM + 1;

    tmpRes(idxM) = gPOCRE(trY,trX,inFamily,inLink,lambda(end),inOption);
    if( bInitW&&(tmpRes(end).nCmp>0) )    % Update W
        fprintf('\b  / ');
        %inOption.W = CalcOptW(tmpRes(end).mu+trX*tmpRes(end).beta,inFamily,inLink);
        eta = bsxfun(@plus,trX*tmpRes(end).beta,tmpRes(end).mu);
        inOption.W = CalcOptW(eta,inFamily,inLink);
        tmpRes(idxM) = gPOCRE(trY,trX,inFamily,inLink,lambda(end),inOption);
    end
    
    %vdn = size(vdX,1);
    %tmpEta = repmat(tmpRes(end).mu,vdn,1)+vdX*tmpRes(end).beta;
    tmpEta = bsxfun(@plus,vdX*tmpRes(end).beta,tmpRes(end).mu);

    if( bMultiN )
        tmpEY = arrayfun(@(i)(ginv(tmpEta(i,:),inLink)),[1:n]','UniformOutput',false);
        tmpEY = cell2mat(tmpEY);
        tmpEY = [1-sum(tmpEY,2) tmpEY];
    else
        tmpEY = ginv(tmpEta,inLink);
    end
    switch( inOption.vc )
        case 'mr'   % miss-classification
            if( strcmp(inFamily,'binomial') )
                yPred = round(tmpEY);
                cvc(idxM) = sum(abs(vdY-yPred))/vdn;
            elseif( bMultiN )
                yPred = bsxfun(@eq,tmpEY,max(tmpEY,[],2));
                cvc(idxM) = cvc(idxM)+sum(sum(abs(tsY-yPred)))/vdn;
            else
                error('  MR cannot be reported for non-binary data!');
            end

        case 'press'    % Default 
            cvc(idxM) = sum(sum((vdY-tmpEY).^2))/vdn;

        case 'loglik'
            if( strcmp(inFamily,'binomial') )
                cvc(idxM) = sum(sum(vdY.*log(tmpEY)+...
                                (1-vdY).*log(1-tmpEY)))/vdn;
            elseif( bMultiN )
                cvc(idxM) = sum(sum(vdY.*log(tmpEY)))/vdn;
            else
                error(' Not Implemented!');
            end
        
        otherwise
            error(['  ' inOption.vc ' is not supported!']);
    end

    nzBeta = sum(abs(tmpRes(end).beta)>eps);
    bSparse = tmpRes(end).bSparse;
    if( lambda(end)>=inOption.lambda0 )
        if( (nzBeta==0)&&bSparse )
            lambda(end+1) = inOption.lambda0-inOption.delta;
            inOption.W = tmpRes(1).W;
        else
            lambda(end+1) = lambda(end)+inOption.delta;
        end
    elseif( (lambda(end)<inOption.lambda0)&&bSparse )
        lambda(end+1) = lambda(end)-inOption.delta;
    else
        break;
    end
end

idxOpt = find(cvc==min(cvc));
idxOpt = idxOpt(1); % In case of multiple minimums
retTP = lambda(idxOpt);
retM = tmpRes(idxOpt);
