function [cvValues,obj] = crossvalidation(obj,numberOfBatches,modelComplexities,ini)
% Evaluates the errors of a v fold cross validation
% (v=numberOfBatches) for the model complexities specified by the user
%
%
% object = crossvalidation(object,numberOfBatches,modelComplexities)
%
%
% Inputs
%
% object:               model object
% numberOfBatches:      integer; in how many batches the whole
%                       dataSet should be devided
% modelComplexities:  	(1 x number of model complexities) vector
%                       containing only integer values; every integer value
%                       determines a model complexity for which the
%                       crossvalidation value should be calculated.
%                       Example: modelComplexities = [3 5] --> The
%                       crossvalidation values for 3 and 5 local models
%                       will be calculated.

% Julian Belz, 16-März-2012
% Institut für Mess- und Regelungstechnik, Universität Siegen, Deutschland
% Institute of Measurement and Control, University of Siegen, Germany
% Copyright (c) 2012 by Julian Belz

if ~exist('numberOfBatches','var') || isempty(numberOfBatches)
    % If there is no numberOfBatches passed, a LOOCV should be performed
    numberOfBatches     = size(obj.input,1);
end

if ~exist('modelComplexities','var') || isempty(modelComplexities)
    % If modelComplexities is empty, perform crossvalidation for the
    % suggested net
    modelComplexities   = obj.suggestedNet;
end

% If the model object contains no data no crossvalidation can be performed
if isempty(obj.input) || isempty(obj.output)
    fprintf('There is no data available to perform a crossvalidation!\n');
    return;
end

if ~exist('ini','var')
    ini = 9596;
end

% information about the size of all data and about the train data size in
% every loop
N           = size(obj.input,1);

% If the demanded number of batches is greater than the number of data
% samples, set the number of batches to the number of data samples
if numberOfBatches-N > eps
    numberOfBatches = N;
    fprintf('Number of demanded batches is greater than the number of data samples.\n Number of batches is set to the number of data samples.\n \n')
end

% Calculate the mean of the number of samples per batch
groupSize   = N - floor(N/numberOfBatches);

% Predefine the crossValidationValues-property of the model-object or
% enlarge this variable
cvValues = zeros(1,max(modelComplexities));

% waitbar
hwait           = waitbar(0,'Please wait...');

% To maintain the training of the original model, the training for the
% crossvalidation calculation is performed with tmpObj
tmpObj          = obj;

% Switch off all termination criterions except for the the number of
% parameters
if isprop(tmpObj,'maxPenaltyDeterioration')
    % Make sure the number of desired models will be reached
    tmpObj.maxPenaltyDeterioration = inf;
end
if isprop(tmpObj,'maxNumberOfLM')
    tmpObj.maxNumberOfLM = inf;
end
if isprop(tmpObj,'minError')
    tmpObj.minError = -inf;
end
if isprop(tmpObj,'maxTrainTime')
    tmpObj.maxTrainTime = inf;
end
if isprop(tmpObj,'maxIterations')
    tmpObj.maxIterations = inf;
end
if isprop(tmpObj,'maxValidationDeterioration')
    tmpObj.maxValidationDeterioration = inf;
end


for jj = [modelComplexities;1:size(modelComplexities,2)]
    
%     % choose current model complexity
%     tmpObj.leafModels = tmpObj.history.leafModelIter{jj(1)};
    
    % set the maximum number of local models for the following
    % cross validation
    tmpObj.maxNumberOfParameters = obj.history.currentNumberOfParameters(jj(1));
    
    % perform cross validation for the current model
    if jj(1) > groupSize
        % If there are more local models allowed than training data
        % points are available, the CV value is set to not a number
        % to avoid problems in the parameter estimation
        cvValues(jj(1)) = NaN;
    else
        Jcv                 = zeros(1,numberOfBatches);
        helpMeCount         = 1;
        output              = Jcv; % Initialising
        modelOutput         = Jcv; % Initialising
        weightings          = Jcv; % Initialising
        
        for ii = 1:numberOfBatches
            % get the samples for training and testing for the current run
            [testSetIdx trainSetIdx]          =  generateCVsamples(obj,numberOfBatches,ii);
            % testSet and trainSet are scaled if the object has been
            % scaled.
            
            % Delete already existing training
            tmpObj.history         = modelHistory;
            tmpObj.leafModels      = [];
            tmpObj.localModels     = [];
            tmpObj.outputModel     = [];
            tmpObj.xRegressor      = [];
            tmpObj.zRegressor      = [];
            if isprop(tmpObj,'phi')
                tmpObj.phi = [];
            end
            if isprop(tmpObj,'MSFValue')
                tmpObj.MSFValue = {[]};
            end
            
            % Train model with new train dataset
            tmpObj.input           = obj.input(trainSetIdx,:);
            tmpObj.output          = obj.output(trainSetIdx,:);
            tmpObj                 = tmpObj.train;
            
            
            % Make sure the correct model complexity is chosen for further steps
            [~,idxMin] = min(abs(tmpObj.history.currentNumberOfParameters - tmpObj.maxNumberOfParameters));
            tmpObj.leafModels      = tmpObj.history.leafModelIter{idxMin};
            
            
            % evaluate the output of the trained model at the test samples
            yTestData               = tmpObj.calculateModelOutput(obj.unscaledInput(testSetIdx,:),obj.unscaledOutput(testSetIdx,:));
            
            
            % evaluation of the squared errors
            outputIndexes           = helpMeCount:(size(yTestData,1)+helpMeCount-1);
            output(1,outputIndexes) = obj.unscaledOutput(testSetIdx,:);
            modelOutput(1,outputIndexes) = yTestData;
            weightings(1,outputIndexes) = obj.dataWeighting(testSetIdx);

            % increasing the variable helpMeCount for the next run of the loop
            helpMeCount             = helpMeCount + size(yTestData,1);
        end
        
        % Calculation of the RMSE
        cvValues(jj(1)) = obj.calcGlobalLossFunction(output', modelOutput', weightings, obj.outputWeighting, obj.lossFunctionGlobal);
    end
    if exist('hwait','var')
        waitbar(jj(2)/size(modelComplexities,2),hwait);
    end
end

% Give only non-zero values back corresponding to the order of the
% 'modelComplexities' vector
cvValues = cvValues(modelComplexities);
obj.numberOfCVgroups = numberOfBatches;
obj.history.kFoldCVlossFunction(modelComplexities) = cvValues;


if exist('hwait','var')
    % hwait does not exist, if the user closed the waitbar window before
    % the calculation ended
    delete(hwait);
end

    function [testSetIdx trainSetIdx] = generateCVsamples(dataSetObject,numberOfGroups,groupSelection,options)
        % Function creates various dataSets
        % [testSet trainSet] = generateTrainTestSet(data, numberOfGroups, groupSelection, options)
        %
        % Output:
        %   testSet:          group size x no. of data columns     test data set
        %   trainSet:         N-group size x no. of data columns   training data set
        %
        % Input:
        %   obj:              dataSet-object
        %   numberOfGroups:   N x q   Number of groups (LOO-cross-validation: numberOfGroups = N)
        %   groupSelection:   1 x 1   Selection of the group that has to be extracted as test set
        %   options:          string  Options for data partitioning
        
        % HILOMOT Nonlinear System Identification Toolbox, Version 1.0
        % Benjamin Hartmann, 27-September-2008
        % Institute of Automatic Control & Mechatronics, Department of Mechanical Engineering, University of Siegen, Germany
        % Copyright (c) 2008 by Benjamin Hartmann
        
        data      = [dataSetObject.input dataSetObject.output];
        
        % Get constants
        Ndata = size(data,1);
        groupsize = floor(Ndata/numberOfGroups);
        
        
        % usualy the groups should be selected randomly
        if nargin < 4
            options = 'randomize';
        end
        
        % Options
        switch options
            case 'randomize'             % Randomize data set
                rng(ini,'twister'); % 9596
                permIdx = randperm(Ndata);
        end
        
        % Perform data partitioning
        if groupSelection == 1
            testSetIdx = permIdx(1:groupsize);
            trainSetIdx = permIdx(groupsize+1:size(data,1));
        elseif groupSelection == numberOfGroups
            testSetIdx = permIdx((groupSelection-1)*groupsize+1 :size(data,1));
            trainSetIdx = permIdx(1:(groupSelection-1)*groupsize);
        else
            testSetIdx = permIdx((groupSelection-1)*groupsize+1 : groupSelection*groupsize);
            trainSetIdx = [permIdx(1:(groupSelection-1)*groupsize), permIdx(groupSelection*groupsize+1:size(data,1))];
        end
    end

end % end crossvalidation