function retRes = SimuStudyCVGPOCRE(inFamily,inLink,modid,nd,n,p,ptype)
%function retRes = SimuStudyCVGPOCRE(inFamily,inLink,modid,nd,n,p,ptype)
%
% Simulation study of gPOCRE for generalized linear models.
%
% Model: g(E[Y]) = mu + X*beta
%
% Input:
%   inFamily: family of the distribution;
%   inLink: link function of the generalized linear model.
%   modid: index of the models;
%   nd: number of datasets;
%   n: number of observation in each dataset;
%   p: number of variables in X;
%   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).
%
% Output:
%   retRes: lambda -- tuning parameter;
%           mr -- misclassification rate;
%           press -- sum of squared prediction error;
%           ntp -- number of true positives (in estimated beta);
%           ntn -- number of true negatives (in estimated beta);
%           nfp -- number of false positives (in estimated beta);
%           nfn -- number of false negatives (in estimated beta);
%           fdr -- false discovery rate (=nfp/(ntp+nfp));
%
% Notes on Family & Link
%   Binomial: logit (default), probit, cloglog
%   Multinomial Baseline Category: bclogit (default)
%   Multinomial Proportional Odds: pologit (default)
%   Multinomial Continuation Ratio: crlogit (default)
%   Poisson: log (default), identity, sqrt
%   Gaussian: identity (default), log, inverse
%   Gamma: inverse (default), identity, log
%   Inverse Gaussian (IGaussian): 1/mu^2 (default), inverse, identity, log
%   Quasi: logit (default), probit, cloglog, identity, inverse, log,
%          1/mu^2, sqrt
%
% Note: quasi family should be able to specify its own power link function
%
% Dabao Zhang    July 1, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(3, 7, nargin));
if( nargin<7 ) ptype = 'ebtz'; end
if( nargin<6 ) p = 1000; end
if( nargin<5 ) n = 100; end
if( nargin<4 ) nd = 100; end

% Initial values
nFolds = 5;    % (training, validation, test) if 1; CV otherwise.
inOption.tol = 1e-6;
inOption.maxit = 50;
inOption.bctype = 1;    % Approximate approach by default
inOption.maxcmp = 50;
inOption.maxvar = n/2;
inOption.ptype = ptype;
%inOption.delta = 0.02;
inOption.delta = 0.1;
inOption.vc = 'mr';
%inOption.vc = 'press';      % Default for validation criterion
retRes.lambda = zeros(nd+3,1);
retRes.mr = zeros(nd+3,1);
retRes.press = zeros(nd+3,1);
retRes.ntp = zeros(nd+3,1);
retRes.ntn = zeros(nd+3,1);
retRes.nfp = zeros(nd+3,1);
retRes.nfn = zeros(nd+3,1);
retRes.fdr = zeros(nd+3,1);

if( isnumeric(inLink) )
    fPrefix = [upper(inFamily(1:3)) num2str(inLink) 'M' num2str(modid) ...
               'N' num2str(n) 'P' num2str(p)];
else
    fPrefix = [upper(inFamily(1:3)) inLink 'M' num2str(modid) ...
               'N' num2str(n) 'P' num2str(p)];
end

for j = 1:nd
    disp([fPrefix ': dataset=' num2str(j)]);
    
    rand('state',j*100);
    randn('state',j*100);
    
    % Generate training data
    tmpYX = SimuSGLMData(inFamily,inLink,modid,2*n+200,p);
    %dataFName = [fPrefix 'DS' num2str(j) '.txt'];
    %YX = [tmpYX.Y tmpYX.X];  save(dataFName,'YX','-ascii'); clear YX;
    
    % Determine the optimal number of components
    if( nFolds==1 ) % training, validation, & test data
        [tmpTP,fRes] = ValidateTP(tmpYX.Y(1:n,:),tmpYX.X(1:n,:),...
                                  tmpYX.Y((n+1):(2*n),:),...
                                  tmpYX.X((n+1):(2*n),:),...
                                  inFamily,inLink,inOption);
        disp(['  Optimal tuning parameter: ' num2str(tmpTP)]);
    else % cross-validation
        inOption.W = [];
        tmpTP = CVGPOCRE(tmpYX.Y(1:n,:),tmpYX.X(1:n,:),inFamily,inLink,...
                         nFolds,inOption);
        disp(['  Optimal tuning parameter: ' num2str(tmpTP)]);
        
        inOption.W = ones(n,1)/n;
        fRes = gPOCRE(tmpYX.Y(1:n,:),tmpYX.X(1:n,:),inFamily,inLink,tmpTP,inOption);
        if( fRes.nCmp>0 )
            % Update W
            fprintf('\b  / ');
            %inOption.W = CalcOptW(fRes.mu+tmpYX.X(1:n,:)*fRes.beta,inFamily,inLink);
            eta = bsxfun(@plus,tmpYX.X(1:n,:)*fRes.beta,fRes.mu);
            inOption.W = CalcOptW(eta,inFamily,inLink);
            fRes = gPOCRE(tmpYX.Y(1:n,:),tmpYX.X(1:n,:),inFamily,inLink,tmpTP,inOption);
        end
    end
    
    % Evaluate the performance
    tsYX.Y = tmpYX.Y((2*n+1):(2*n+200),:);
    tsYX.X = tmpYX.X((2*n+1):(2*n+200),:);
    pfmRes = EvalPerf4GLM(inFamily,inLink,tmpYX.beta,fRes,tsYX);
    
    % Save the results
    retRes.lambda(j) = tmpTP;
    if( ~isempty(pfmRes.mr) ) retRes.mr(j) = pfmRes.mr; end
    if( ~isempty(pfmRes.press) ) retRes.press(j) = pfmRes.press; end
    retRes.ntp(j) = pfmRes.ntp;
    retRes.ntn(j) = pfmRes.ntn;
    retRes.nfp(j) = pfmRes.nfp;
    retRes.nfn(j) = pfmRes.nfn;
    retRes.fdr(j) = pfmRes.fdr;
end

% Calculate summary statistics: mean, median, standard deviation
retRes.lambda(nd+1) = mean(retRes.lambda(1:nd),1);
retRes.lambda(nd+2) = median(retRes.lambda(1:nd),1);
retRes.lambda(nd+3) = std(retRes.lambda(1:nd));
retRes.mr(nd+1) = mean(retRes.mr(1:nd),1);
retRes.mr(nd+2) = median(retRes.mr(1:nd),1);
retRes.mr(nd+3) = std(retRes.mr(1:nd));
retRes.press(nd+1) = mean(retRes.press(1:nd),1);
retRes.press(nd+2) = median(retRes.press(1:nd),1);
retRes.press(nd+3) = std(retRes.press(1:nd));
retRes.ntp(nd+1) = mean(retRes.ntp(1:nd),1);
retRes.ntp(nd+2) = median(retRes.ntp(1:nd),1);
retRes.ntp(nd+3) = std(retRes.ntp(1:nd));
retRes.ntn(nd+1) = mean(retRes.ntn(1:nd),1);
retRes.ntn(nd+2) = median(retRes.ntn(1:nd),1);
retRes.ntn(nd+3) = std(retRes.ntn(1:nd));
retRes.nfp(nd+1) = mean(retRes.nfp(1:nd),1);
retRes.nfp(nd+2) = median(retRes.nfp(1:nd),1);
retRes.nfp(nd+3) = std(retRes.nfp(1:nd));
retRes.nfn(nd+1) = mean(retRes.nfn(1:nd),1);
retRes.nfn(nd+2) = median(retRes.nfn(1:nd),1);
retRes.nfn(nd+3) = std(retRes.nfn(1:nd));
retRes.fdr(nd+1) = mean(retRes.fdr(1:nd),1);
retRes.fdr(nd+2) = median(retRes.fdr(1:nd),1);
retRes.fdr(nd+3) = std(retRes.fdr(1:nd));

% Save the results
resFName = [fPrefix 'T' upper(ptype) 'CVRes.csv'];
csvwrite(['Lambda-' resFName],retRes.lambda);
csvwrite(['MR-' resFName],retRes.mr);
csvwrite(['PRESS-' resFName],retRes.press);
csvwrite(['NTP-' resFName],retRes.ntp);
csvwrite(['NTN-' resFName],retRes.ntn);
csvwrite(['NFP-' resFName],retRes.nfp);
csvwrite(['NFN-' resFName],retRes.nfn);
csvwrite(['FDR-' resFName],retRes.fdr);
      
disp(['    MR:    median=' num2str(retRes.mr(nd+2))]);
disp(['    MR:        sd=' num2str(retRes.mr(nd+3))]);
disp(['    PRESS: median=' num2str(retRes.press(nd+2))]);
disp(['    PRESS:     sd=' num2str(retRes.press(nd+3))]);
