function pred = DNNTEST(optnet,optlrtheta, cv_data, cv_label)
% data : mxp training data matrix (no need to normalize), m is # of samples and p is # of dims
% label: training labels \in {0,1}, can be binary (class) or real-valued (reg)
% cv_data: xp dev data matrix, used for early stopping
% cv_label: corresponding dev data labels

% neural net model consists of 'optnet','optlrtheta','para'. So just save
% those.DNNTraina.m
%5-3 YJ change to test only DNNTEST.M

layer1=optnet(1);
hiddL1n=size(layer1.W,1);
layer2=optnet(2);
hiddL2n=size(layer2.W,1);

%% params setting
para.hidstruct = [hiddL1n hiddL2n]; % num of hidden layers and units

%para.isRBM = true; % RBM pretraining? For small dataset, this will significantly improve performance
%para.isRBM = false;
para.rbmEpoch = 60; 

para.sgdBatch = 256*1; % # of samples in a minibatch
para.sgdMaxEpoch = 50; % num of minibatch gradient descent for backprop
para.lambda_L2 = 1e-5;

% para.thetaType = 'avg';
para.thetaType = 'hitfa'; % early stopping criterion: accuracy, hit-fa, average

para.isMVA = true; % MVA preprocessing?
%para.isMVA = false;
para.MVAOrder = 2; % order of ARMA filtering in MVA preprocessing

para.costFunc = 'xentropy'; % cross entropy
% para.costFunc = 'maxHF'; % hit-fa maximization
% para.costFunc = 'HFAcc'; % cross+entropy + hit-fa
% para.costFunc = 'consFA'; % hit-fa max with constraint on fa
%% normalization
if para.isMVA
   [cv_data, para.tr_mu, para.tr_std] = MVAProcess(cv_data, para.MVAOrder, 0);
    %[rowTrData para.tr_mu para.tr_std] = MVAProcess(data, para.MVAOrder, 0); 
    %cv_data = MVAProcess_testing(cv_data, para.MVAOrder, 0, para.tr_mu, para.tr_std);%shouldmake sure;jy
else
    [cv_data, para.tr_mu, para.tr_std] = mean_var_norm(data);
    %[rowTrData, para.tr_mu, para.tr_std] = mean_var_norm(data);
    %cv_data = mean_var_norm_testing(cv_data, para.tr_mu, para.tr_std);
end
%{
%% intialization
netStruct = [size(rowTrData,2) para.hidstruct]
if para.isRBM
%     preNet = trainStackRBM(rowTrData, netStruct, para.rbmEpoch, 512);
    preNet = trainStackRBM(rowTrData, netStruct, para.rbmEpoch, 256);
else
    preNet = randinitnet(netStruct,1);
end
    
preopttheta = 0.01*randn(netStruct(end)+1,1);
%}
%% CV
cv_data(isnan(cv_data)) = 0;
cv_label(cv_label == -1) = 0;
%% SGD training
%{
maxEpoch = para.sgdMaxEpoch;
epsilonw = linspace(1,0.001,maxEpoch); % linearly decrease learning rates

initialmomentum  = 0.5;
finalmomentum    = 0.9;

label(label== -1) = 0 ;
nettheta = unRolling(preNet);
alltheta = [nettheta; preopttheta];

batchid = genBatchFixNum(size(data,1),para.sgdBatch);
nBatch = size(batchid,2);

para
Winc = 0;
icount = 0;
avgtheta = 0;
cvrecord = [];
tic
for epoch = 1: maxEpoch
   seq = randperm(size(data,1));      
   permData = rowTrData(seq,:);
   permLabel = label(seq,:);
   
   costsum = 0;
   for bid = 1: nBatch-1 % need to fix it
       icount = icount + 1;
       numcases = batchid(2,bid) - batchid(1,bid) + 1;
       batchData = permData(batchid(1,bid):batchid(2,bid),:);
       batchLabel = permLabel(batchid(1,bid):batchid(2,bid));
       
       switch para.costFunc
           case 'xentropy'
               [cost grad] = fineTuneLRCost(alltheta, netStruct, batchLabel, batchData, para.lambda_L2);
           case 'maxHF'
               [cost grad] = fineTuneHFLRCost(alltheta, netStruct, batchLabel, batchData, para.lambda_L2);
           case 'HFAcc'
               [cost grad] = fineTuneHFAccLRCost(alltheta, netStruct, batchLabel, batchData, para.lambda_L2, 0.5);               
           case 'consFA'
               [cost grad] = fineTuneHFLRCost_consFA(alltheta, netStruct, batchLabel, batchData, para.lambda_L2, 0.05,0.8);                                             
           case 'alter'
%                if mod(epoch,2)
               if epoch <= 8
                   [cost grad] = fineTuneLRCost(alltheta, netStruct, batchLabel, batchData, para.lambda_L2);
               else
                   [cost grad] = fineTuneHFLRCost(alltheta, netStruct, batchLabel, batchData, para.lambda_L2);
               end
           otherwise
               error('no cost func.')
       end
      
       if epoch>maxEpoch/2
         momentum=finalmomentum;
       else
         momentum=initialmomentum;
       end
      
       Winc = momentum*Winc + epsilonw(epoch)*(grad);       
       alltheta = alltheta - Winc;                     
       costsum = costsum + cost;
   end   
   disp(['Sum of cost at epoch ', int2str(epoch), ' is ', num2str(costsum)]);   
   
   avgtheta = avgtheta + alltheta;
   
   % cross-validation 
   [hitfa_rate acc_rate] = checkCVRate(alltheta, cv_data, cv_label, netStruct);
   disp(['CV HIT-FA Rate: ' num2str(hitfa_rate) ' Acc Rate: ' num2str(acc_rate)])
   cvrecord.hitfa_rate(epoch) = hitfa_rate;
   cvrecord.acc_rate(epoch) = acc_rate;
   cvrecord.X(epoch,:) = alltheta';   
end
SGDTime = toc

switch para.thetaType
    case 'avg'
        optalltheta = avgtheta / maxEpoch;
        [hitfa_rate acc_rate hit fa cvpred] = checkCVRate(optalltheta, cv_data, cv_label, netStruct);
        disp(['Avergared SGD --> CV HIT-FA Rate: ' num2str(hitfa_rate) ' Acc Rate: ' num2str(acc_rate)])        
    case 'hitfa'
        [mv mi] = max(cvrecord.hitfa_rate);
        disp(['Early Stopping --> best HIT-FA rate ' num2str(mv),' at epoch ',int2str(mi)]);
        optalltheta = cvrecord.X(mi,:)';
        [hitfa_rate acc_rate hit fa cvpred] = checkCVRate(optalltheta, cv_data, cv_label, netStruct);
        disp(['HIT: ' num2str(hit) ' FA: ' num2str(fa)])   
        outputhit=hit;
        outputfa=fa;
        outputacc=acc_rate;
        outputpred=cvpred;%label out
    case 'acc'
        [mv mi] = max(cvrecord.acc_rate);
        disp(['Early Stopping --> best Acc rate ' num2str(mv),' at epoch ',int2str(mi)]);
        optalltheta = cvrecord.X(mi,:)';                
        [hitfa_rate acc_rate hit fa] = checkCVRate(optalltheta, cv_data, cv_label, netStruct);
        disp(['HIT: ' num2str(hit) ' FA: ' num2str(fa)])
end

[optnet optlrtheta] = Rolling(optalltheta, netStruct);
%}
%%test 
% change 'rowTrData' to your test data to get prediction of your test data (remeber to normalize like we did before)
[lastact] = passThruNet(optnet, cv_data); 
[pred, prob]= lrRegressionPred(optlrtheta,lastact); % prob is posterior prob, sometimes more useful than simple binary labels
%acc = mean(pred == label);
pred = double(pred);

