function obj = train(obj)
% TRAIN is the central training algorithm of the HiLoMoT toolbox. A given data set is is modeled
% with local polynomial models. The interpolation between the local models is accomplished with
% validity functions that are hierarchically constructed with sigmoidal splitting functions. The
% splitting with HILOMOT is performed as long as the maximum number of local models is not
% achieved or the training error doesn't go below a user defined limit, respectively.
%
%
% obj = train(obj)
%
%
% INPUT:
%
%   obj             object   HILOMOT object containing all relevant properties and methods.
%
%
% OUTPUT:
%
%   obj             object   HILOMOT object containing all relevant properties and methods.
%
%
% HiLoMoT - Nonlinear System Identification Toolbox
% Benjamin Hartmann, 04-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles

% Check regressors for rule premises (z) and rule consequents (x)
if isempty(obj.xRegressor) || isempty(obj.zRegressor)
    obj.xRegressor = data2xRegressor(obj,obj.input,obj.output);
    obj.zRegressor = data2zRegressor(obj,obj.input,obj.output);
end

% Set timer for training
tStart = tic;

% Estimate first net
if isempty(obj.localModels) % Estimate a global linear model
    obj = estimateFirstLM(obj); % one linear model
else % Estimate a nonlinear model based on the given structure
    warning('hilomot:train','Training based on a given model is not implemented so far. Training aborted!')
    return
end

% Write history
obj.history = writeHistory(obj.history,obj,tStart);

% Display loss function value for the first net
if obj.history.displayMode
    fprintf('\nInitial net has %d local model(s): J = %.3f.\n', ...
        numel(obj.localModels), obj.history.globalLossFunction(end));
end

% Tree construction algorithm
while obj.checkTerminationCriterions % Go on training until one termination criterion is reached
    
    while any(obj.idxAllowedLM & obj.leafModels) % While-loop over all allowed local models
        
        % Initialize the flag, if a permitted initial split is found
        flagPermittedInitialSplitFound = false;
        
        % Initialize the flag, if a permitted split is found
        flagPermittedSplitFound = false;
        
        % Find worst performing LM
        worstLM = findWorstLM(obj);
        
        % Check for split of worst performing LLM
        if obj.history.displayMode
            fprintf('\n\n%d. Iteration. Number of local models = %d. Checking for split of model %d ...', ...
                obj.history.iteration, sum(obj.leafModels), worstLM);
        end
        
        % Initialize the error for the best starting point of the optimization
        globalLossFunctionBestIni = inf;
        
        % Initialize the error for the best model, in every iteration the
        % training error must decrease
        globalLossFunctionBest = obj.history.globalLossFunction(end);
        
        % Perform axes-orthogonal splitting
        for splitDimension = 1:size(obj.zRegressor,2) % Over all dimensions of the z-regressor
            [objSplitOrtho, forbiddenSplit] = LMSplit(obj, worstLM, splitDimension);
            
            % If the split is forbidden, try the next dimension or split ratio
            if forbiddenSplit
                if obj.history.displayMode
                    fprintf('\n   Testing split in dimension %d:         Forbidden!', splitDimension);
                end
                continue % Continue with the next split
            end
            
            if obj.history.displayMode
                fprintf('\n   Testing split in dimension %d:         J = %.3f.', splitDimension, objSplitOrtho.history.globalLossFunction(end));
            end
            
            % If better than the currently best axes-orthogonal split
            if globalLossFunctionBestIni - objSplitOrtho.history.globalLossFunction(end) > eps 
                objBestIni = objSplitOrtho;
                globalLossFunctionBestIni = objSplitOrtho.history.globalLossFunction(end);
                flagPermittedInitialSplitFound = true; % At least one permitted othogonal split is found
            end
        end
        
        % Perform axes-oblique splitting
        if obj.oblique
            
            % Test previous split for possibly better initialization
            if sum(obj.leafModels)>1 && size(obj.zRegressor,2) > 1
                splitDimension = 'lastSplit'; % Over all dimensions of the z-regressor
                [objSplitLast, forbiddenSplit] = LMSplit(obj, worstLM, splitDimension);
                
                if forbiddenSplit
                    if obj.history.displayMode
                        fprintf('\n   Testing previous split:               Forbidden!');
                    end
                else
                    if obj.history.displayMode
                        fprintf('\n   Testing previous split:               J = %.3f.', objSplitLast.history.globalLossFunction(end));
                    end
                    % If better than the best axes-orthogonal split
                    if objSplitLast.history.globalLossFunction(end) < globalLossFunctionBestIni
                        objBestIni = objSplitLast;
                        globalLossFunctionBestIni = objSplitLast.history.globalLossFunction(end);
                        flagPermittedInitialSplitFound = true; % At least the last split direction leads to a permitted split
                    end
                end
            end
            
            % Check if a permitted initial split exist; without initialization no nonlinear optimization
            if flagPermittedInitialSplitFound
                
                % Set best initial model as currently best model, if it is better the model of the last iteration
                if globalLossFunctionBestIni < globalLossFunctionBest % Check if the initial model is better
                    objBest = objBestIni;
                    globalLossFunctionBest = globalLossFunctionBestIni;
                    flagPermittedSplitFound = true; % A permitted split is found, that decrease the training error
                end
                
                % Split optimization
                [objSplitOblique, forbiddenSplit] = LMSplitOpt(obj, objBest.localModels(end-1).splittingParameter, worstLM);
                
                if forbiddenSplit
                    if obj.history.displayMode
                        fprintf('\n   Axes-oblique splitting:               Forbidden!');
                    end
                else
                    if obj.history.displayMode
                        fprintf('\n   Axes-oblique splitting:               J = %.3f.', objSplitOblique.history.globalLossFunction(end));
                    end
                    
                    % If better than the model of the last iteration and the best initial model
                    if objSplitOblique.history.globalLossFunction(end) < globalLossFunctionBest  % Check if the optimized model is better
                        objBest = objSplitOblique;
                        flagPermittedSplitFound = true; % A split was found that decrease the training error
                    end
                end
            end
            
        elseif flagPermittedInitialSplitFound % Check if one permitted axes-orthogonal split exist
            % Check if the orthogonal split leads to animprovement of the performance
            if objBestIni.history.globalLossFunction(end) < globalLossFunctionBest
                objBest = objBestIni;
                flagPermittedSplitFound = true; % A permitted split is found, that decrease the training error
            end
        end
        
        % Check if a permitted split of the current worstLM can be done
        if flagPermittedSplitFound
            % Break while-loop; no other local models need to be tested
            break
        else
            % No permitted split found; go one investigating another local model and lock the worstLM
            obj.idxAllowedLM(worstLM) = false;
            if obj.history.displayMode
                fprintf('\n\n%d. Iteration. Model %d can not be splitted. Try another model...', ...
                    obj.history.iteration, worstLM);
            end
        end
        
    end
    
    % Update model object, if a permitted split is found
    if flagPermittedSplitFound
        
        % Perform best split
        obj = objBest;
        
        % Make history entries
        obj.history = writeHistory(obj.history,obj,tStart);
        
        if obj.history.displayMode
            fprintf('\n-> SPLITTING RESULT:                     J = %.3f.', obj.history.globalLossFunction(end));
        end
    end
    
end

% Perform reliability check of new model
obj.numberOfLMReliable = sum(obj.idxAllowedLM & obj.leafModels);

% Final display
if obj.history.displayMode
    fprintf('\n\nFinal net has %d local models and %d parameters: J = %.3f.', ...
        sum(obj.leafModels), round(obj.history.currentNumberOfParameters(end)), obj.history.globalLossFunction(end));
    fprintf('\n\nNet %d with %d LMs and %d parameters is suggested as the model with the best complexity trade-off.\n\n', ...
        obj.suggestedNet, sum(obj.history.leafModelIter{obj.suggestedNet}), round(obj.history.currentNumberOfParameters(obj.suggestedNet)));
end
end






