function retRes = SimuStudyGPOCREPath(inFamily,inLink,modid,nd,n,p,ptype,inTheta)
%function retRes = SimuStudyGPOCREPath(inFamily,inLink,modid,nd,n,p,ptype,inTheta)
%
% Simulation study of tuning parameter selection criteria when applying
% gPOCRE for generalizeg 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);
%   inTheta: signal strength (=1 by default);
%
% 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
%   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    May 21, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(3, 8, nargin));
if( nargin<8 ) inTheta = 1; end
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
%inOption.tol = 1e-8;
inOption.tol = 1e-6;
inOption.maxit = 50;
inOption.bctype = 1;    % Approximate approach by default
inOption.maxcmp = 50;
%inOption.maxvar = n/2;
inOption.maxvar = n/log(n); % Follow SIS by Fan (????)
inOption.ptype = ptype;
%inOption.gamma = 10;   % =3.7 by default
%inOption.delta = 0.02; % =0.1 by default)
%inOption.lambda0 = 1.2;    % =1 by default
retRes.lambda = zeros(nd+3,13);
retRes.mr = zeros(nd+3,13);
retRes.press = zeros(nd+3,13);
retRes.ntp = zeros(nd+3,13);
retRes.ntn = zeros(nd+3,13);
retRes.nfp = zeros(nd+3,13);
retRes.nfn = zeros(nd+3,13);
retRes.fdr = zeros(nd+3,13);

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

% 2 ~ AIC; 3 ~ AICc; 0:0.1:1 ~ (E)BIC
msc = [2 3 0:0.1:1];
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,n+200,p,inTheta);
    dataFName = [fPrefix 'DS' num2str(j) '.txt'];
    YX = [tmpYX.Y tmpYX.X];  save(dataFName,'YX','-ascii'); clear YX;
    tmpRes = gPOCREPath(tmpYX.Y(1:n,:),tmpYX.X(1:n,:),inFamily,inLink,inOption);
    
    % Test Data
    tsYX.Y = tmpYX.Y((n+1):end,:);
    tsYX.X = tmpYX.X((n+1):end,:);
    
    % MSC
    for k=1:13
        fRes = SelectModel(tmpRes,msc(k));
        pfmRes = EvalPerf4GLM(inFamily,inLink,tmpYX.beta,fRes,tsYX);
        retRes.lambda(j,k) = fRes.lambda;
        if( ~isempty(pfmRes.mr) ) retRes.mr(j,k) = pfmRes.mr; end
        if( ~isempty(pfmRes.press) ) retRes.press(j,k) = pfmRes.press; end
        retRes.ntp(j,k) = pfmRes.ntp;
        retRes.ntn(j,k) = pfmRes.ntn;
        retRes.nfp(j,k) = pfmRes.nfp;
        retRes.nfn(j,k) = pfmRes.nfn;
        retRes.fdr(j,k) = pfmRes.fdr;
    end
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' num2str(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,:))]);
