function [PreBootStrap, PostBootStrap] =  GetPerformanceINRIA(options, savename);
%function [PreBootStrap, PostBootStrap] =  GetPerformanceINRIA(options, savename);
%
%computes performance on the INRIA Dataset
%assumes standards on the data formatting and filename conventions
%specifically, PosTrainFN = sprintf('PosTrainData_%s.mat', FeatureName);
%NegTrainFN = sprintf('PosTrainData_%s_Part%.3d.mat', FeatureName, PartNumber);
%
%Assume we always use all the positive training and testing data
%what remains to be decided is what part of the negative data will be used for
%  training, bootstraping, and testing

HelpTextInfo = 'GetPerformanceINRIA';
AppRoot = '/cbcl/scratch01/bileschi/HoughTForm';
seedrand;

if(nargin < 1), options= [], end;
d.Verbose = 1;
d.DataRoot = '/cbcl/scratch03/bileschi/Images/INRIAPerson/MatrixFormat';
d.cFeatureString{1} = 'GS';
d.cVarName{1} = 'X';
d.nNegTrainingSamples = 4000;
d.nPosTrainingSamples = inf;
d.nNegTestingSamples = 16000;
d.bDoBootstrapping = 1;
d.nBootstrapSamples = 50000;
d.nMaxBootstapSamples = 4000;
d.sPARAMS.Nrounds = 100;
d.sPARAMS.R = 0;
d.sPARAMS.T = 0;
d.sPARAMS.C = 1;
d.sPARAMS.KERNEL = 0;
d.nMaxFeatures = 4000;
d.ClassifierName = 'gentleBoost';
d.featuresubset = [];
d.RandSeed = floor(rand*99999);
d.bNormalizeDataSamples = 1;
d.caNormalizationRegions = [];
d.NormalizationRelativeStrengths = [1];
d.bAbsVal = 0;
d.nNegTestPartsAvailable = 50;
d.nNegTrainPartsAvailable = 50;
options = ResolveMissingOptions(options,d);

fprintf('Building Training Data... ');
[yTrainFull, XTrainFull,options] = GetDataX(options, 'Train');
[XTrainFull, FeatureSubsample] = SubSampleFeatures(XTrainFull, options);
[XTrainFull,NormalizationParams] = InternalNormalize(XTrainFull, options);
fprintf('%d x %d\n', size(XTrainFull,1), size(XTrainFull,2));
fprintf('Training Classifier\n');
yTrainFull = (yTrainFull(:))';
cmd = sprintf('Model_1 = CLS%s(XTrainFull,yTrainFull'',options.sPARAMS);',options.ClassifierName);
eval(cmd);
L = length(yTrainFull);

fprintf('Testing Classifier... ');
[yel_1, ye_1,yTestFull] = INRIATestClassifier(options, FeatureSubsample, Model_1,NormalizationParams);


ROC_1 = ROCetc2(ye_1',yTestFull);
PreBootStrap.ROC = ROC_1;
PreBootStrap.Model = Model_1;
PreBootStrap.ye = ye_1;
PreBootStrap.ytest = yTestFull;
PostBootStrap = [];
if(options.bDoBootstrapping)
  fprintf('Building Bootstrap Data... ');
  nBS = 0;
  for i = 1:ceil(options.nBootstrapSamples / 1000)
    [yBSFull, XBSFull] = GetDataX(options, 'Bootstrap', i);
    XBSFull = SubSampleFeatures(XBSFull, options, FeatureSubsample);
    XBSFull = InternalNormalize(XBSFull, options,NormalizationParams);
    fprintf('%d x %d\n', size(XTrainFull,1), size(XTrainFull,2));
    fprintf('Testing Bootstrap Data group %d of %d... ', i, ceil(options.nBootstrapSamples / 1000));
    cmd = sprintf('[yel, ye] = CLS%sC(XBSFull, Model_1);',options.ClassifierName);
    eval(cmd);
    iNewExmpls = find(yel > 0);
    nBS = nBS + length(iNewExmpls);
    fprintf(' %d new samples selected\n', length(iNewExmpls));
    XTrainFull = [XTrainFull, XBSFull(:,iNewExmpls)]; 
    yTrainFull = [yTrainFull, yBSFull(iNewExmpls)];
    if(nBS > d.nMaxBootstapSamples)
      break;
    end 
  end
  clear XBSFull;
  fprintf('Training Bootstrapped Classifier\n');
  cmd = sprintf('Model_2 = CLS%s(XTrainFull,yTrainFull'',options.sPARAMS);',options.ClassifierName);
  eval(cmd);
  % Model_2 = CLSgentleBoost(XTrainFull,yTrainFull',options.sPARAMS);
  fprintf('Testing Bootstrapped Classifier\n');
  [yel_2, ye_2,yTestFull] = INRIATestClassifier(options, FeatureSubsample, Model_2,NormalizationParams);
  ROC_2 = ROCetc2(ye_2',yTestFull);
  PostBootStrap.ROC = ROC_2;
  PostBootStrap.Model = Model_2;
  PostBootStrap.ye = ye_2;
  PostBootStrap.ytest = yTestFull;
end
if(exist('savename'))
  save(savename,'HelpTextInfo','PreBootStrap','PostBootStrap','options');
end


function [yfull, NewDataX,options] = GetDataX(options,Epoch, FileMandate)
NewDataX = [];
yfull = [];
nCellParts = length(options.cVarName);
nParts.Train = ceil(options.nNegTrainingSamples / 1000); % 1000 samples per file
nParts.Test = ceil(options.nNegTestingSamples / 1000); % 1000 samples per file
% nParts.Bootstrap = ceil(options.nNegBootstrapSamples / 1000); % 1000 samples per file
for iFeature = 1:nCellParts
  switch lower(Epoch)
    case {'train'}
      [cyPos{iFeature}, cXPos{iFeature}] = GetDataXHelper(options, 1,  ...
        sprintf('PosTrainData_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature});
      [cyNeg{iFeature}, cXNeg{iFeature},p] = GetDataXHelper(options, nParts.Train , ...
        sprintf('NegTrainData_Part*_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature});
      options.UsedTrainParts = p;
      if(size(cXPos{iFeature},2) > options.nPosTrainingSamples)
        seedrand(options.RandSeed);
        p = randperm(size(cXPos{iFeature},2));
        p = p(1:options.nPosTrainingSamples);
        cXPos{iFeature} = cXPos{iFeature}(:,p);
	cyPos{iFeature} = cyPos{iFeature}(p);
      end
    case {'bootstrap'}
      cyPos{iFeature} = [];
      cXPos{iFeature} = [];
      [cyNeg{iFeature}, cXNeg{iFeature},p] = GetDataXHelper(options, 1,  ...
        sprintf('NegTrainData_Part*_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature}, FileMandate);
    case {'test'}
      if(exist('FileMandate'))
        cyPos{iFeature} = [];
	cXPos{iFeature} = [];
	[cyNeg{iFeature}, cXNeg{iFeature},p] = GetDataXHelper(options, 1, ...
          sprintf('NegTestData_Part*_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature}, FileMandate);
        options.UsedTestParts = FileMandate;
      else
        [cyPos{iFeature}, cXPos{iFeature}] = GetDataXHelper(options, 1, ...
          sprintf('PosTestData_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature});
        [cyNeg{iFeature}, cXNeg{iFeature},p] = GetDataXHelper(options, nParts.Test, ...
          sprintf('NegTestData_Part*_%s.mat',options.cFeatureString{iFeature}), options.cVarName{iFeature});
        options.UsedTestParts = p;
      end
  end
  NewDataX = [NewDataX;[cXPos{iFeature},cXNeg{iFeature}]];
end % (loop over multiple features)
if(options.bAbsVal)
  NewDataX = abs(NewDataX);
end
yfull = [cyPos{1},cyNeg{1}];
usedNegParts = p;

function [y,X,p] = GetDataXHelper(options,nParts,Template,Varname,FileMandate);
d = dir(fullfile(options.DataRoot,Template));
nPartsAvail = length(d);
if(nParts > nPartsAvail)
  fprintf('WARNING: REQUESTED TOO MANY PARTS, USING ALL PARTS\n');
  nParts = nPartsAvail;
end
if(nPartsAvail == 0 )
  error('ERROR: No Data Available.  Check Filename\n%s\n',fullfile(options.DataRoot,Template));
end
seedrand(options.RandSeed);  
p = randperm(nPartsAvail);
p = p(1:nParts);
X = [];
y = [];
if(nargin > 4)
  p = intersect(FileMandate,1:nPartsAvail);
end
for iFl = p;
  D = load(fullfile(options.DataRoot, d(iFl).name));%--> y, (Varname) 
  yt = getfield(D,'y');
  Xt = getfield(D,Varname);
  if(isfield(options,'featuresubset'))
    if(not(isempty(options.featuresubset)))
      Xt = Xt(options.featuresubset,:);
    end
  end
  X = [X,Xt];
  y = [y,yt];
end

function [XTrainFull, SelectedFeatures] = SubSampleFeatures(XTrainFull, options, SelectedFeatures);
if(nargin < 3)
  SelectedFeatures = 1:size(XTrainFull,1);
  if(size(XTrainFull,1) > options.nMaxFeatures), 
    seedrand(options.RandSeed);
    p = randperm(size(XTrainFull,1));
    XTrainFull = XTrainFull(p(1:options.nMaxFeatures),:);
    SelectedFeatures = p(1:options.nMaxFeatures);
    fprintf('TOO MANY FEATURES: Randomly sampling %d features\n', options.nMaxFeatures);
  end
else
  XTrainFull = XTrainFull(SelectedFeatures,:);
end  

function [yel, ye, ytru] = INRIATestClassifier(options, FeatureSubsample, Model,NormalizationParams)
% In order to not have to hold all the data in memory at once, we load the negatives in parts, 
% and store only the classification score.
TempnNegTestingSamples = options.nNegTestingSamples;
nParts.Test = ceil(options.nNegTestingSamples / 1000); % 1000 samples per file
% d = dir(fullfile(options.DataRoot,sprintf('NegTestData_Part*_%s.mat',options.cFeatureString{iFeature})));
% nPartsAvail = length(d);
nPartsAvail = options.nNegTestPartsAvailable;
seedrand(options.RandSeed);
p = randperm(nPartsAvail);
p = p(1:nParts.Test);
options.nNegTestingSamples = 0;
[ytru, XTestPos] = GetDataX(options,'Test');
XTestPos = SubSampleFeatures(XTestPos, options, FeatureSubsample);
XTestPos = InternalNormalize(XTestPos,options,NormalizationParams);
cmd = sprintf('[yel, ye] = CLS%sC(XTestPos,Model);',options.ClassifierName);
eval(cmd);
n = 1;
for i = p
  fprintf('Testing Classifier, part %d of %d\n',n,length(p));
  [yt, XtNeg] = GetDataX(options,'Test', i);
  XtNeg = SubSampleFeatures(XtNeg, options, FeatureSubsample);
  XtNeg = InternalNormalize(XtNeg,options,NormalizationParams);
  cmd = sprintf('[yel_t, ye_t] = CLS%sC(XtNeg,Model);',options.ClassifierName);
  eval(cmd);
  yel = [yel(:); yel_t(:)];
  ye = [(ye(:))', (ye_t(:))'];
  ytru = [(ytru(:))', (yt(:))'];
  n = n+1;
end
fprintf('\n');
options.UsedTestParts = p;
options.nNegTestingSamples = TempnNegTestingSamples;

function [X,NormalizationParams] = InternalNormalize(X,options,NormalizationParams)
if(nargin < 3)
  NormalizationParams = [];
end
if(options.bNormalizeDataSamples > 0)
  if(isempty(options.caNormalizationRegions))
    caNR{1} = 1:size(X,1);
  else
    caNR = options.caNormalizationRegions;
  end
  [X,NormalizationParams] = normalizeByIndex(X, caNR, options.bNormalizeDataSamples,NormalizationParams,options);
end
     
