function retRes = CompareGPOCRE2GNET(inFamily,modid,nd,n,p,ptype)
%function retRes = CompareGPOCRE2GNET(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, 9, 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

addpath('./glmnet');

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
    
% 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.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);

% Set up the store matrix for comparision method
ResComp.mr = zeros(nd+3,1);
ResComp.ntp = zeros(nd+3,1);
ResComp.ntn = zeros(nd+3,1);
ResComp.nfp = zeros(nd+3,1);
ResComp.nfn = zeros(nd+3,1);
ResComp.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) ...
        'ND' num2str(nd) 'N' num2str(n) 'P' num2str(p) 'ptype' ptype];
end

% Set up directory for storing the summary statistics
mkdir(fPrefix);

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);
	%save the generated dataset.
    dataFName = [fPrefix 'DS' num2str(j) '.txt'];
    YX = [tmpYX.Y tmpYX.X];  save(['./' fPrefix '/' 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);
    
    % Extract the training dataset
    tranYX.Y = tmpYX.Y(1:n,:);
    tranYX.X = tmpYX.X(1:n,:);
    % Evaluate the performance of glmnet
    switch inFamily
        case 'multibc'
            retResComp=EvalPerf4GNET_BC(tranYX,tmpYX.beta,tsYX,nFolds);
        case 'multipo'
            retResComp=EvalPerf4GNET_PO(tranYX,tmpYX.beta,tsYX,nFolds );
        case 'multicr'
            retResComp=EvalPerf4GNET_CR(tranYX,tmpYX.beta,tsYX,nFolds);
        otherwise
            error('  SimuStudyCVPOCRE: The family is not supported in multinomail model ');
    end
    % Save the results for the comparision method
    if( ~isempty(retResComp.mr) ) ResComp.mr(j) = retResComp.mr; end
    ResComp.ntp(j) = retResComp.ntp;
    ResComp.ntn(j) = retResComp.ntn;
    ResComp.nfp(j) = retResComp.nfp;
    ResComp.nfn(j) = retResComp.nfn;
    ResComp.fdr(j) = retResComp.fdr;
    
    % Save the results for POCRE
    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 for POCRE
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 for pocre
resFName = [fPrefix 'T' upper(ptype) 'CV_POCRE.csv'];
csvwrite(['./' fPrefix '/Lambda-' resFName],retRes.lambda);
csvwrite(['./' fPrefix '/MR-' resFName],retRes.mr);
csvwrite(['./' fPrefix '/PRESS-' resFName],retRes.press);
csvwrite(['./' fPrefix '/NTP-' resFName],retRes.ntp);
csvwrite(['./' fPrefix '/NTN-' resFName],retRes.ntn);
csvwrite(['./' fPrefix '/NFP-' resFName],retRes.nfp);
csvwrite(['./' fPrefix '/NFN-' resFName],retRes.nfn);
csvwrite(['./' fPrefix '/FDR-' resFName],retRes.fdr);

% Calculate summary statistics: mean, median, standard deviation for
% for comparision method

ResComp.mr(nd+1) = mean(ResComp.mr(1:nd),1);
ResComp.mr(nd+2) = median(ResComp.mr(1:nd),1);
ResComp.mr(nd+3) = std(ResComp.mr(1:nd));
ResComp.ntp(nd+1) = mean(ResComp.ntp(1:nd),1);
ResComp.ntp(nd+2) = median(ResComp.ntp(1:nd),1);
ResComp.ntp(nd+3) = std(ResComp.ntp(1:nd));
ResComp.ntn(nd+1) = mean(ResComp.ntn(1:nd),1);
ResComp.ntn(nd+2) = median(ResComp.ntn(1:nd),1);
ResComp.ntn(nd+3) = std(ResComp.ntn(1:nd));
ResComp.nfp(nd+1) = mean(ResComp.nfp(1:nd),1);
ResComp.nfp(nd+2) = median(ResComp.nfp(1:nd),1);
ResComp.nfp(nd+3) = std(ResComp.nfp(1:nd));
ResComp.nfn(nd+1) = mean(ResComp.nfn(1:nd),1);
ResComp.nfn(nd+2) = median(ResComp.nfn(1:nd),1);
ResComp.nfn(nd+3) = std(ResComp.nfn(1:nd));
ResComp.fdr(nd+1) = mean(ResComp.fdr(1:nd),1);
ResComp.fdr(nd+2) = median(ResComp.fdr(1:nd),1);
ResComp.fdr(nd+3) = std(ResComp.fdr(1:nd));

% Save the results for glmnet
resFName = [fPrefix 'T' upper(ptype) 'CV_GLMNET.csv'];
csvwrite(['./' fPrefix '/MR-' resFName],ResComp.mr);
csvwrite(['./' fPrefix '/NTP-' resFName],ResComp.ntp);
csvwrite(['./' fPrefix '/NTN-' resFName],ResComp.ntn);
csvwrite(['./' fPrefix '/NFP-' resFName],ResComp.nfp);
csvwrite(['./' fPrefix '/NFN-' resFName],ResComp.nfn);
csvwrite(['./' fPrefix '/FDR-' resFName],ResComp.fdr);

% Save summary statistics
methods = char('glmnet','pocre');
smr = char([num2str(ResComp.mr(nd+1)) '(' num2str(ResComp.mr(nd+2)) ')'],...
    [num2str(retRes.mr(nd+1)) '(' num2str(retRes.mr(nd+2)) ')']);
sfdr = char([num2str(ResComp.fdr(nd+1)) '(' num2str(ResComp.fdr(nd+2)) ')'],...
    [num2str(retRes.fdr(nd+1)) '(' num2str(retRes.fdr(nd+2)) ')']);
sntp = char([num2str(ResComp.ntp(nd+1)) '(' num2str(ResComp.ntp(nd+2)) ')'],...
    [num2str(retRes.ntp(nd+1)) '(' num2str(retRes.ntp(nd+2)) ')']);
sntn = char([num2str(ResComp.ntn(nd+1)) '(' num2str(ResComp.ntn(nd+2)) ')'],...
    [num2str(retRes.ntn(nd+1)) '(' num2str(retRes.ntn(nd+2)) ')']);
snfp = char([num2str(ResComp.nfp(nd+1)) '(' num2str(ResComp.nfp(nd+2)) ')'],...
    [num2str(retRes.nfp(nd+1)) '(' num2str(retRes.nfp(nd+2)) ')']);
snfn = char([num2str(ResComp.nfn(nd+1)) '(' num2str(ResComp.nfn(nd+2)) ')'],...
    [num2str(retRes.nfn(nd+1)) '(' num2str(retRes.nfn(nd+2)) ')']);
SummaryTable = table(smr,sfdr,snfp,snfn,sntp,sntn);
sresFName = [fPrefix 'T' upper(ptype) 'CV_GLMNET_POCRE_Summary.csv'];
writetable(SummaryTable,['./' fPrefix  '/' sresFName]);
