function retTP = CVGPOCRE(inY,inX,inFamily,inLink,nFolds,inOption)
%function retTP = CVGPOCRE(inY,inX,inFamily,inLink,nFolds,inOption)
%
% Optimize the number of components using cross-validation for Penalized
% Orthogonal-Component Regression in fitting generalized linear models
% (gPOCRE) with high dimensional data. Please refer details in
%
% Lin, Y., Zhang, M., and Zhang, D. (2011) Penalized orthogonal-component
% regression for high-dimensional generalized linear models. Purdue
% Statistics Technical Report ???.
%
% Model: g(E[Y}) = mu + X*beta;
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inFamily: the distribution family, such as 'normal', 'binomial',
%             'poisson', 'gamma', 'inverse gaussian', 'quasi' (not
%             supported now);
%   inLink: the link function g(.), such as 'identity' (default for
%           'normal'), 'log' (default for 'poisson'), 'logit' (default for
%           'binomial'), 'probit','comploglog','reciprocal','loglog'
%           (default for 'gamma'),p (numeric, default for 'inverse
%           gaussian'), user-defined link function (not supported now);
%   nFolds: number of folds for the cross-validation (=10 by default);
%   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 (gamma=3.7 by default),
%                'mcp' ~ MCP (gamma=3.7 by default);
%       gamma -- parameter for SCAD or MCP (=3.7 by default);
%       delta -- pace for choosing value of the tuning parameter;
%       lambda0 -- initial value of the tuning parameter (1 by default);
%       vc -- cross-validation criterion: 'press' (default),'mr','loglik';
%
% Output:
%   retTP: optimal value of the tuning parameter.
%
% 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(3, 6, nargin));

if( nargin<6 ) inOption = []; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; 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<5 ) nFolds = 10; end

if( (nargin<4)||isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        otherwise
            error('  CVGPOCRE: The family is not supported...');
    end
end

[n,p] = size(inX);    % X with n objects p variables;
k = size(inY,2);    % Y with n objects k variables;
if( n~=size(inY,1) )
    error('  CVGPOCRE: inX and inY should have the same number of rows!');
end

% Only implemented for single response!
if( k>1 )
    error('  CVGPOCRE: Only implemented for single response!');
end

if( nFolds<2 )
    warning('  CVGPOCRE: The number of folds is reset to 2 ...');
    nFolds = 2;
end

% Randomly split the data into nCVFold folds
rndIdx = randperm(n);
sgmIdx = ceil([1:n]/round(n/nFolds));
if( any(sgmIdx>nFolds) )
    tmpIdx = find(sgmIdx>nFolds);
    sgmIdx(tmpIdx) = 1:length(tmpIdx);
    sgmIdx = sort(sgmIdx);
end

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

    % Go through all folds for the same lambda:
    cvc(idxM) = 0;
    nzBeta = 0;
    for j=1:nFolds
        disp(['    fold = ' num2str(j)]);

        % training dataset
        trX = inX(rndIdx(sgmIdx~=j),:);
        trY = inY(rndIdx(sgmIdx~=j),:);

        % test dataset
        tsX = inX(rndIdx(sgmIdx==j),:);
        tsY = inY(rndIdx(sgmIdx==j),:);
        [rowtsY,coltsY]=size(tsY);
        
        inOption.W = W(rndIdx(sgmIdx~=j));
        tmpRes = gPOCRE(trY,trX,inFamily,inLink,lambda(end),inOption);
        if( bInitW&&(tmpRes.nCmp>0) ) % Update W if necessary
            fprintf('\b  / ');
            inOption.W = CalcOptW(tmpRes.mu+trX*tmpRes.beta,inFamily,inLink);
            tmpRes = gPOCRE(trY,trX,inFamily,inLink,lambda(end),inOption);
        end
        
        tmpEta = repmat(tmpRes.mu,rowtsY,1)+tsX*tmpRes.beta;
        tmpEY = ginv(tmpEta,inLink);
        %tmpEY = ginv(bsxfun(@plus,tsX*tmpRes.beta,tmpRes.mu),inLink);
        switch( inOption.vc )
            case 'mr'   % miss-classification
                if( strcmp(inFamily,'binomial') )
                    yPred = round(tmpEY);
                    cvc(idxM) = cvc(idxM)+sum(abs(tsY-yPred))/rowtsY;
                else
                    error('  MR cannot be reported for non-binary data!');
                end
                
            case 'press'    % Default 
                cvc(idxM) = cvc(idxM)+sum(sum((tsY-tmpEY).^2))/rowtsY;
                
            case 'loglik'
                cvc(idxM) = sum(sum(tsY.*log(tmpEY)+...
                                    (1-tsY).*log(1-tmpEY)))/rowtsY;
            
            otherwise
                error(['  ' inOption.vc ' is not supported!']);
        end
        
        if( j==1 )
            nzBeta = sum(abs(tmpRes.beta)>eps);
            bSparse = tmpRes.bSparse;
        else
            nzBeta = min(nzBeta, sum(abs(tmpRes.beta)>eps));
            bSparse = max(bSparse, tmpRes.bSparse);
        end
    end
    
    if( lambda(end)>=inOption.lambda0 )
        if( (nzBeta==0)&&bSparse )
            lambda(end+1) = inOption.lambda0-inOption.delta;
        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

retTP = lambda(cvc==min(cvc));
retTP = retTP(1); % in case of multiple minimums
