function [PreBootStrap, PostBootStrap] =  GetPerformanceINRIA_RealBootstrap(options, savename);
%function [PreBootStrap, PostBootStrap] =  GetPerformanceINRIA_RealBootstrap(options, savename);
%

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

if(nargin < 1), options= [], end;
d.Verbose = 1;
d.PositiveDataRoot = '/cbcl/scratch03/bileschi/Images/INRIAPerson/MatrixFormat';
d.NegativeDataRoot = '/cbcl/scratch03/bileschi/Images/INRIAPerson/WindowingFormat';
d.cFeatureString{1} = 'GS32x16';
d.cVarName{1} = 'X';
d.nNegTrainingSamples = 500;
d.nMaxBootstrapSamples = 5000;
d.sPARAMS.Nrounds = 100;
d.sPARAMS.R = 0;
d.sPARAMS.T = 0;
d.sPARAMS.C = 1;
d.sPARAMS.KERNEL = 0;
d.nMaxFeatures = 2500;
d.ClassifierName = 'libsvm';
d.featuresubset = [];
d.RandSeed = floor(rand*99999);
d.bNormalizeDataSamples = 1;
d.caNormalizationRegions = [];
d.featureSpatialyx = [32 16];
d.ExamplewiseNormalization = 'L2';
d.bNormalizeFeatures = 0;
d.FeaturewiseNormalization = 'outnormalize';
d.FeatureNormalizationParams = [];
d.bUseLocalNeighborhoodSuppression = 0;

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... ');
% INRIATestClassifier needs to be edited to use the windowing framework instead.
[yel_1, ye_1,yTestFull] = INRIATestClassifier(options, FeatureSubsample, Model_1,NormalizationParams);

ROC_1 = ROCetc2(ye_1(:),yTestFull(:),1);
PreBootStrap.ROC = ROC_1;
PreBootStrap.Model = Model_1;
PreBootStrap.ye = ye_1;
PreBootStrap.ytest = yTestFull;
PostBootStrap = [];
if(exist('savename'))
  save(sprintf('%s_PreBootOnly.mat',savename),'HelpTextInfo','PreBootStrap','options');
end
% Model one Complete, Begin Bootstrapping stage.

[XBootstrap, yBootstrap] = INRIAGetBootstrapExamples(options, FeatureSubsample, Model_1, NormalizationParams);
XTrainFull = [XTrainFull, XBootstrap]; 
yTrainFull = [yTrainFull(:)', yBootstrap(:)']; 
fprintf('Training Bootstrapped Classifier\n');
cmd = sprintf('Model_2 = CLS%s(XTrainFull,yTrainFull'',options.sPARAMS);',options.ClassifierName);
eval(cmd);
fprintf('Testing Bootstrapped Classifier\n');
[yel_2, ye_2,yTestFull] = INRIATestClassifier(options, FeatureSubsample, Model_2,NormalizationParams);
ROC_2 = ROCetc2(ye_2',yTestFull,1);
PostBootStrap.ROC = ROC_2;
PostBootStrap.Model = Model_2;
PostBootStrap.ye = ye_2;
PostBootStrap.ytest = yTestFull;
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 = inf; % 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;
    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} = [];
        % [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)
yfull = [cyPos{1},cyNeg{1}];

function [y,X,p] = GetDataXHelper(options,nParts,Template,Varname,FileMandate);
d = dir(fullfile(options.PositiveDataRoot,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.PositiveDataRoot,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.PositiveDataRoot, 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 [XBootstrap, yBootstrap] = INRIAGetBootstrapExamples(options, FeatureSubsample, Model, NormalizationParams);
seedrand(options.RandSeed);
NegTestRoot = '/cbcl/scratch03/bileschi/Images/INRIAPerson/WindowingFormat/Train';
FeatString = options.cFeatureString{1};
if(strcmpi(FeatString,'GS32x16'))
  FeatString = 'GS';
end
d = dir(fullfile(NegTestRoot,FeatString,'*.mat'));
nfiles = length(d);
nWindows = 0;
nf = min(options.nMaxFeatures, length(FeatureSubsample))
XBootstrap = zeros(nf, options.nMaxBootstrapSamples);
yeB = repmat(-inf,[1,size(XBootstrap,2)]);
for imgIdx = 1:nfiles
  fprintf('Boot Boot Boot DATA Boostrap File %d of %d \n',imgIdx, nfiles);
  load(fullfile(NegTestRoot, FeatString,d(imgIdx).name)); %--> cImgPyr;
  [X,y] = INRIAGetBootstrapSamplesOnePyr(cImgPyr, Model,options,FeatureSubsample,NormalizationParams);
  [XBootstrap, yeB] = PickHardestExamples(XBootstrap,yeB, X, y, options);
end      
yBootstrap = -1 + 0 * yeB;

function  [XBootstrap, yeB] = PickHardestExamples(XBootstrap,yeB, X, y, options);
My = [yeB(:);y(:)];
MX = [XBootstrap, X];
[s,si] = sort(-My);
si = si(1:min(length(si),options.nMaxBootstrapSamples));
XBootstrap = MX(:,si);
yeB = My(si);  


function [yel, ye, ytru] = INRIATestClassifier(options, FeatureSubsample, Model,NormalizationParams)
% Postive testing examples are pre extracted.
% Negative testing occurs as a windowing process over a set of pre built image pyramids.
% Each negative window of size options.featureSpatialyx is tested.
seedrand(options.RandSeed);
[ytru, XTestPos] = GetDataX(options,'Test'); % --> only positive data
XTestPos = SubSampleFeatures(XTestPos, options, FeatureSubsample);
XTestPos = InternalNormalize(XTestPos, options,NormalizationParams);
cmd = sprintf('[yel_pos, ye_pos] = CLS%sC(XTestPos,Model);',options.ClassifierName);
eval(cmd);
NegTestRoot = '/cbcl/scratch03/bileschi/Images/INRIAPerson/WindowingFormat/Test';
FeatString = options.cFeatureString{1};
if(strcmpi(FeatString,'GS32x16'))
  FeatString = 'GS';
end
d = dir(fullfile(NegTestRoot,FeatString,'*.mat'));
nfiles = length(d);
nWindows = 0;
nFP = zeros(length(ye_pos),1);
for imgIdx = 1:nfiles
  fprintf('TEST TEST TEST DATA Test File %d of %d \n',imgIdx, nfiles);
  load(fullfile(NegTestRoot, FeatString,d(imgIdx).name)); %--> cImgPyr;
  [nWindows_t, nFP_t] = INRIAGetModelPerformanceOnOnePyr(cImgPyr, Model,ye_pos,options,FeatureSubsample,NormalizationParams);
  nWindows = nWindows + nWindows_t;
  nFP = nFP + nFP_t;
end      
[yel, ye, ytru] = ConstructPlausibleDetections(nWindows, nFP, ye_pos);

function [X,y] = INRIAGetBootstrapSamplesOnePyr(Pyr, Model,options,FeatureSubsample,NormalizationParams);
nPyrLevels = length(Pyr);
options.ReturnBootstrapExamples = 1;
options.MaxBootstrapExamples = 5000;
options.FeatureSubsample = FeatureSubsample;
options.NormalizationParams = NormalizationParams;
X = [];
y = [];
for i = 1:nPyrLevels
  [yEmpl,yEmp,XBootStrap,yEmpBootStrap] = WindowedObjectDetection2(Model, options.featureSpatialyx, Pyr{i},options);
  [X, y] = PickHardestExamples(XBootStrap,yEmpBootStrap, X, y, options);
end
 

function [nWindows, nFP] = INRIAGetModelPerformanceOnOnePyr(Pyr,Model,ye,options,FeatureSubsample,NormalizationParams);
%for every threshold, this returns the number of windows scanned in the pyramid, and the number of
%detections within those windows. 
nPyrLevels = length(Pyr);
options.ReturnBootstrapExamples = 0;
options.MaxBootstrapExamples = 5000;
s = sort(ye);
l = length(s);
nWindows = 0;
nFP = zeros(size(s));
options.FeatureSubsample = FeatureSubsample;
options.NormalizationParams = NormalizationParams;
for i = 1:nPyrLevels
  [yEmpl,yEmp,XBootStrap,yEmpBootStrap] = WindowedObjectDetection2(Model, options.featureSpatialyx, Pyr{i},options);
  if(options.bUseLocalNeighborhoodSuppression)
     [yEmp, yEmpl] = ApplyLocalNeighborhoodSuppression(yEmp, yEmpl);
  end
  nWindows = nWindows + length(yEmp(:));
  for j = 1:length(s)
    nFP(j) = nFP(j) + length(find(yEmp(:) >= s(j)));
  end
end

function [yel, ye, ytru] = ConstructPlausibleDetections(nWindows, nFP, ye_pos);
%build yel, ye, and ytru such that they match the nWindows, nFP and ye_pos given.
%this will allow us to still use ROCetc2 rather than writing new ROC and DET code.
ye = zeros(nWindows+length(ye_pos),1);
ye_pos = sort(ye_pos);
surrogateNegScore = (ye_pos(1:(end-1)) + ye_pos(2:end)) / 2;
surrogateNegScore = [ye_pos(1) - 1, surrogateNegScore', ye_pos(end)+1];
n_of_each_neg(1) = nWindows - nFP(1);
for i = 2:length(nFP)
   n_of_each_neg(i) = - nFP(i) + nFP(i-1);
end
n_of_each_neg(length(nFP) + 1) = nFP(length(nFP));
n = 0
for i = 1:length(surrogateNegScore)
  ye(n + (1:n_of_each_neg(i))) = surrogateNegScore(i);
  n = n+n_of_each_neg(i);
end   
ye((n+1):end) = ye_pos;
ytru = [-ones(1,nWindows),ones(1,length(ye_pos))];
yel = ye > 0;

