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;
%       Dabao Zhang, Mar. 5, 2016 -- Correct the case idxM>inOption.maxit.
%
% 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 '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('  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

% Special for multinomial models
bMultiN = 0;
if( sum(strcmp(inFamily,{'multibc','multipo','multicr'})) )
    bMultiN = 1;
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);
            eta = bsxfun(@plus,trX*tmpRes.beta,tmpRes.mu);
            inOption.W  = CalcOptW(eta,inFamily,inLink);
            tmpRes = gPOCRE(trY,trX,inFamily,inLink,lambda(end),inOption);
        end
        
        %tmpEta = repmat(tmpRes.mu,rowtsY,1)+tsX*tmpRes.beta;
        tmpEta = bsxfun(@plus,tsX*tmpRes.beta,tmpRes.mu);
        
        if( bMultiN )
            tmpEY = arrayfun(@(i)(ginv(tmpEta(i,:),inLink)),[1:rowtsY]',...
                             '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) = cvc(idxM)+sum(abs(tsY-yPred))/rowtsY;
                elseif( bMultiN )
                    yPred = bsxfun(@eq,tmpEY,max(tmpEY,[],2));
                    cvc(idxM) = cvc(idxM)+sum(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;
                if( strcmp(inFamily,'binomial') )
                    cvc(idxM) = sum(sum(tsY.*log(tmpEY)+...
                                    (1-tsY).*log(1-tmpEY)))/rowtsY;
                elseif( bMultiN )
                    cvc(idxM) = sum(sum(tsY.*log(tmpEY)))/rowtsY;
                else
                    error(' Not Implemented!');
                end
                                
            otherwise
                error(['  ' inOption.vc ' is not supported!']);
        end
        
        if( j==1 )
            nzBeta = sum(sum(abs(tmpRes.beta)>eps));
            bSparse = tmpRes.bSparse;
        else
            nzBeta = min(nzBeta, sum(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
lambda  = lambda(1:length(cvc));

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

% Q: Should we handle the last non-sparse lambda differently?
%if( (retTP==lambda(end))&&(~bSparse) )
%    retTP = retTP(end-1);
%end
