function retPV = SimuStudyRBootsGPOCRE(inFamily,inLink,modid,nd,n,p,ptype)
%function retPV = SimuStudyRBootsGPOCRE(inFamily,inLink,modid,nd,n,p,ptype)
%
% Simulation study of the generalized Penalized Orthogonal-Components
% REgression (gPOCRE) with p-values calculated from the bootstrap method.
%
% 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, 1~4;
%   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 (default),
%            'l1' ~ L1,
%            'scad' ~ SCAD (with gamma=3.7),
%            'mcp' ~ MCP (with gamma=3.7).
%
% Output:
%   retPV: bootstrap p-values.
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(2, 8, nargin));

if( nargin<8 ) inTheta = 1; end
if( nargin<7 ) ptype = 'ebtz'; end
if( nargin<6 ) p = 1000; end
if( nargin<5 ) n = 200; end
if( nargin<4 ) nd = 100; end
if( nargin<3 ) mid = 1; end

p = max(1000,p);
nBoots = 100;
fPrefix = ['M' num2str(modid) 'N' num2str(n) 'P' num2str(p)];

% Initial values
inOption.ptype = ptype;
inOption.delta = 0.01;

retPV = ones(p,nd);
%retMNZBeta = zeros(p,nd);
%retSDNZBeta = zeros(p,nd);
%retNNZBeta = zeros(p,nd);

for j = 1:nd
    disp([fPrefix ': dataset=' num2str(j)]);
    
    rand('state',j*100);
    randn('state',j*100);
    
    % Generate training data
    simYX = SimuSGLMData(inFamily,inLink,modid,n,p,inTheta);
    %dataFName = [fPrefix 'DS' num2str(j) '.txt'];
    %YX = [tmpYX.Y tmpYX.X];  save(dataFName,'YX','-ascii'); clear YX;
    
    % Screen first
    inOption.maxvar = n;
    inOption.maxcmp = 5;
    [tmpX,tmpXIdx] = gPOCREScreen(simYX.Y,simYX.X,inFamily,inLink,[],inOption);

    % Fit the model
    inOption.maxvar = n/2;
    inOption.maxcmp = 10;
    tmpRes = gPOCREPath(simYX.Y,tmpX,inFamily,inLink,inOption);
    tmpRes = SelectModel(tmpRes);
    
    % Bootstrap to calculate p-values
    tmpRes = RBootsGPOCRE(simYX.Y,tmpX,tmpRes,inFamily,inLink,nBoots,inOption);
    retPV(tmpXIdx,j) = tmpRes.pvalue;
    %retMNZBeta(tmpXIdx,j) = tmpRes.mnzbeta;
    %retSDNZBeta(tmpXIdx,j) = tmpRes.sdnzbeta;
    %retNNZBeta(tmpXIdx,j) = tmpRes.nnzbeta;
end

% Save the results
%resFName = [fPrefix 'T' upper(ptype) 'Res.mat'];
%save(['PV-' resFName],'retPV');
dlmwrite(['PV-' fPrefix 'T' upper(ptype) 'Res.txt'],retPV);
      
disp([fPrefix ': It is done!']);
