function [mean_err,std_err] = nestedCrossValidation(X, y, f_learn, f_class, d, e, lambda)
% perform triple cross validation and report the avg. error
% X is a n,m matrix with n instances with m features
% Y is a n,1 vector with n labels for the n instances
% f_learn is the function handle for the learn function
% f_class is the function handle for the learn function
% d is the number of outer folds
% e is the number of inner folds
% lambda is a vector of used regularization parameter





% number of used training examples per fold
m_cv = floor(size(X,1)/d);

prediction_error = zeros(d,1);
% shuffle the train data
rIdx = randperm(size(X,1));
X = X(rIdx,:);
y = y(rIdx);
for i = 1:d
    % index of test examples
    p_te = (i-1)*m_cv+1:i*m_cv;
    % index of train examples
	p_tr = [1:(i-1)*m_cv i*m_cv+1:size(X,1)];
    
    Xtrain = X(p_tr,:);
    Ytrain = y(p_tr);
    Xtest = X(p_te,:);
    Ytest = y(p_te);
    % number used training examples per fold (for find optimal lambda)
    n_cv = floor(size(Xtrain,1)/e);
    lambda_acc = zeros(length(lambda),e);
    for j = 1:e
        % index of test examples
        p_test = (j-1)*n_cv+1:j*n_cv;
        % index of train examples
        p_train = [1:(j-1)*n_cv j*n_cv+1:size(Xtrain,1)];
        XtrainEval = Xtrain(p_train,:);
        YtrainEval = Ytrain(p_train);
        XtestEval = Xtrain(p_test,:);
        YtestEval = Ytrain(p_test);      
        for k = 1:length(lambda)
			model = f_learn(XtrainEval, YtrainEval, lambda(k));
            Yeval = f_class(model,XtestEval);
            lambda_acc(k,j) = sum(Yeval == YtestEval)/ length(YtestEval);
        end
    end
    avg_accuracy = sum(lambda_acc,2);
    [maxVal,maxIdx] = max(avg_accuracy);    
    model = f_learn(Xtrain,Ytrain,lambda(maxIdx));
    Yeval = f_class(model,Xtest);
    prediction_error(i) = 1 - (sum(Ytest == Yeval) / length(Yeval));
    disp(['iteration: ' num2str(i) ' opt C: ' num2str(lambda(maxIdx)) ' error: ' num2str(prediction_error(i))]);
end

mean_err = mean(prediction_error);
std_err = std(prediction_error) / sqrt(numel(prediction_error));
disp(['average error: ' num2str(err) ' std-err: ' num2str(std_err)]);
	

end