function [e_crit, prediction, accuracyme] = conformity(TrainX,TrainY,TestX,ratio,gamma,C,TestY)

% Runs a yf(x) conformity measure on SVM with Gaussian kernel 
%   -- requires LibSVM
% 
% Use:
%   [e_crit, prediction] = conformity(TrainX,TrainY,TestX,ratio,gamma,C)
%   [e_crit, prediction, accuracyme] = ... 
%                    conformity(TrainX,TrainY,TestX,ratio,gamma,C,TestY)
%
% Input:
%   TrainX - Training data [samples x features]
%   TrainY - Training data labels [samples x 1]
%   TestX  - Testing data [samples x features]
%   ratio  - the % of training data to be used as calibration data [0<r<1]
%   gamma  - the Gaussian width parameter, inputted as a vector
%   C      - the SVM penalty term, inputted as a vector
%   TestY  - Testing data labels [optional] [samples x 1]
%
% Output:
%   e_crit     - the probability of error (confidence) on the predictions
%   prediction - the prediction values
%   accuracyme - if TestY is inputted then mean accuracy is computed otherwise 0
%
% Written by David R. Hardoon, University College London
%  01/02/09 D.Hardoon@cs.ucl.ac.uk
%
%  Cite: "A Nonconformity Approach to Model Selection for SVMs" by David R.
%  Hardoon, Zakria Hussain and John Shawe-Taylor. University College London, 
%  Dept. of Computer Science 2009
%  http://www.davidroihardoon.com/academic/Publications_files/nonconformity_TR.pdf
%
% Updated:
%   01/04/09 by Zakria Hussain
%   30/06/09 by David R. Hardoon - better error handling
%

%% Error handling
if (nargin < 6 | nargin > 7)
   e_crit = 0; prediction = 0; accuracyme = 0; 
   help conformity
   disp('!!!!! - Incorrect number of inputs - !!!!!');
   return
elseif (nargin == 6)
    TestY = [];
end

if (ratio <= 0 || ratio >= 1)
    e_crit = 0; prediction = 0; accuracyme = 0;
    help conformity
    disp('!!!!! - Can not have Calibration set smaller, larger');
    disp('or equal to training data. Choose a 0 < ratio < 1 - !!!!!');
    return
end

%% Parameter initilisation
num = length(gamma)*length(C);
L = length(TrainY);
accuracyst = 0;
accuracyme = 0;

% Randomising the order of the training data -- so we choose a random
% calibration set
rp = randperm(L);
TrainY = TrainY(rp);
TrainX = TrainX(rp,:);

% Finding out how many positive and negative data we have
npos = sum(TrainY==1);
nneg = sum(TrainY==-1);

% Seperating the positive and negative from the main variable
Xpos = TrainX(TrainY==1,:);
Xneg = TrainX(TrainY==-1,:);

% Selecting the number of positive + negative to be used for calibration
numcal = round((size(TrainX,1)*ratio)/2);

% Setting calibration set
Yval = [ones(numcal,1); -ones(numcal,1)];
Xval = [Xpos(1:numcal,:); Xneg(1:numcal,:)];
nval = length(Yval);

% Updating training data by removing calibration set
Ytrain = [ones(npos-numcal,1); -ones(nneg-numcal,1)];
Xtrain = [Xpos(numcal+1:end,:); Xneg(1+numcal:end,:)];

% Some more variable initialitaion
ntest = size(TestX,1);
count = 1;
val_pos = [];
val_neg = [];
DV_TEST = [];
DV_VAL = [];
val_pos1 = [];
val_neg1 = [];
tmpLabel = ones(size(TestX,1),1);  %% dummy variable.

%% Computing all SVM models
for C_count=1:length(C)
    for gamma_count=1:length(gamma)
        % Training the SVM with C and gamma
        libsvm_options = sprintf('-g %f -c %f', gamma(gamma_count) ,C(C_count));    
        models(C_count,gamma_count).model = svmtrain(Ytrain,Xtrain,libsvm_options);

        % Computing the prediction yf(x) on the validation data
        [PreLabels_val, accuracy_val, decision_values_val] = svmpredict(Yval,Xval, models(C_count,gamma_count).model);
        decision_values_val = PreLabels_val.*abs(decision_values_val);

        % Computing the prediction f(x) on the test data given the current model
        [PreLabels, accuracy, decision_values] = svmpredict(tmpLabel,TestX, models(C_count,gamma_count).model);            
        decision_values = PreLabels.*abs(decision_values);

        % Computing the conformity measure as in paper A(z_j) <= A(z_i) for y = +1
        val_pos(:,count) = (sum(repmat(Yval,1,ntest).*repmat(decision_values_val,1,ntest) <= repmat(decision_values',nval,1))+1)/(nval+1);
        
        % Computing the conformity measure as in paper A(z_j) <= A(z_i) for y = -1
        val_neg(:,count) = (sum(repmat(Yval,1,ntest).*repmat(decision_values_val,1,ntest) <= repmat(-decision_values',nval,1))+1)/(nval+1); 

        % test point decision values
        DV_TEST(:,count) = decision_values;
        
        % decision values for validation points 
        DV_VAL(:,count) = decision_values_val;     

        % Count of number of models
        count = count+1;
    end
end

%% Selection of models with lowest probability and associated label predictions

% validation errors
err = sum(DV_VAL.*repmat(Yval,1,num)<0)/nval;

% sorting the errors (we only want to take the top best models)
[a,b] = sort(err);

% only using the top 20% of best models
L1 = round(num*0.5);
indices = b(1:L1);  
r = randperm(L1)';

for i=1:L1  % number of functions ot use 
    for k=1:size(DV_TEST,1)
        val_pos(k,i) = (sum(DV_VAL(:,indices(r(i))).*Yval <= DV_TEST(k,indices(r(i))))+1)/(nval+1);
        val_neg(k,i) = (sum(DV_VAL(:,indices(r(i))).*Yval <= -DV_TEST(k,indices(r(i))))+1)/(nval+1);     
    end
    model(i,1) = indices(r(i,1));
end

% initialising probabilities
e_crit = zeros(ntest,1);

% random prediction for ties
prediction = sign(randn(ntest,1));  
indx = find(min(val_neg,[],2) == min(val_pos,[],2));
pval = min(val_neg,[],2);
e_crit(indx) = pval(indx);

% predict positive
indx = find(min(val_neg,[],2)< min(val_pos,[],2));
pval = min(val_neg,[],2); 
prediction(indx) = 1;
e_crit(indx) = pval(indx);

% predict negative
indx = find(min(val_pos,[],2)< min(val_neg,[],2));
pval = min(val_pos,[],2); 
prediction(indx) = -1;
e_crit(indx) = pval(indx);

% compute accuracy
if ~isempty(TestY)
    accuracyme = mean(prediction == TestY);
end