function [obj, forbiddenSplit] = LMSplitEstimate(obj, worstLM, splitDimension, splitRatio)
% LMSPLITESTIMATE splits one local model into two halves and
% estimates all relevant properties of the new local model objects.
% Afterwards it evaluates the local and global performance and updates the
% LOLIMOT object.
%
% obj = LMSplitEstimate(obj, worstLM, splitDimension, splitRatio)
%
% OUTPUT:
%   obj:             object     LOLIMOT object containing all relevant net
%                               and data set information and variables.
%
% INPUT:
%   obj:             object     LOLIMOT object containing all relevant net
%                               and data set information and variables.
%   worstLLM:        (1 x 1)    Index of the worst performing LLM which is
%                               splitted.
%   splitDimension:  (1 x 1)    Splitting dimension.
%   splitRatio:      (1 x 1)    Splitting ratio.
%
%
% LoLiMoT - Nonlinear System Identification Toolbox
% Torsten Fischer, 08-December-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


% Set flag
forbiddenSplit = false;

% Get some constants
numberOfLM = numel(obj.leafModels);

% Update active models and allocation of gaussian
obj.leafModels(worstLM) = false;
obj.leafModels = [obj.leafModels  true(1,2)];

% Update the allowed LM
obj.idxAllowedLM = [obj.idxAllowedLM true(1,2)];

% Update net history
obj.history.splitLM = [obj.history.splitLM; worstLM];
obj.history.splitDimension = [obj.history.splitDimension; splitDimension];
obj.history.splitRatio = [obj.history.splitRatio; splitRatio];

% Perform split
[lowerLeftCorner1, upperRightCorner1, lowerLeftCorner2, upperRightCorner2]  = ...
    obj.LMSplit(obj.localModels(worstLM).lowerLeftCorner, obj.localModels(worstLM).upperRightCorner, splitDimension, splitRatio);

% Evaluate centers and standard deviations for the new local models and add the to the lolimot object
obj.localModels = [obj.localModels, gaussianOrthoLocalModel(lowerLeftCorner1,upperRightCorner1,obj.smoothness),...
    gaussianOrthoLocalModel(lowerLeftCorner2,upperRightCorner2,obj.smoothness)];

% Add the training boundaries of the z-regression
obj.localModels(numberOfLM+1).zUpperBound = obj.localModels(worstLM).zUpperBound;
obj.localModels(numberOfLM+1).zLowerBound = obj.localModels(worstLM).zLowerBound;
obj.localModels(numberOfLM+2).zUpperBound = obj.localModels(worstLM).zUpperBound;
obj.localModels(numberOfLM+2).zLowerBound = obj.localModels(worstLM).zLowerBound;

% Evaluate new validity function values
if ~isempty(strfind(version,'R2010a')) 
    % workaround for R2010a
    obj.MSFValue = [obj.MSFValue, obj.localModels(end-1).calculateMSF(obj.zRegressor), obj.localModels(end).calculateMSF(obj.zRegressor)];
else
    obj.MSFValue = [obj.MSFValue, arrayfun(@(loc) loc.calculateMSF(obj.zRegressor),obj.localModels(end-1:end),'UniformOutput',false)];
end

% Get all active gaussians and normalize them
validityFunctionValue = obj.calculateVFV(obj.MSFValue(obj.leafModels));

% Indicies of the active models
idxLeafModels = find(obj.leafModels);

% Check for sufficient amount of data samples
if obj.pointsPerLMFactor * size(obj.xRegressor,2) - sum(validityFunctionValue{end-1}) > eps
    forbiddenSplit = true;
    return
end
if obj.pointsPerLMFactor * size(obj.xRegressor,2) - sum(validityFunctionValue{end}) > eps
    forbiddenSplit = true;
    return
end

% Estimate parameters of the LM
obj.localModels(numberOfLM+1).parameter = estimateParametersLocal(obj,obj.xRegressor, obj.output, validityFunctionValue{end-1}, [], obj.dataWeighting);
obj.localModels(numberOfLM+2).parameter = estimateParametersLocal(obj,obj.xRegressor, obj.output, validityFunctionValue{end}, [], obj.dataWeighting);

% Calculate model output
obj.outputModel = calculateModelOutput(obj);

% Calculate the local loss function values for all leaf models
localLossFunctionValue = calcLocalLossFunction(obj, obj.unscaledOutput, obj.outputModel , cell2mat(validityFunctionValue), obj.dataWeighting, obj.outputWeighting);
for k = 1:numel(idxLeafModels)
    % Allocate the right local loss function value to each local model
    obj.localModels(idxLeafModels(k)).localLossFunctionValue = localLossFunctionValue(k);
end

% Calculate the global loss function value
obj.history.globalLossFunction(end+1) = calcGlobalLossFunction(obj, obj.unscaledOutput, obj.outputModel, obj.dataWeighting, obj.outputWeighting);

end