function [obj, forbiddenSplit] = LMSplit(obj, worstLM, splitDimension)
% LMSPLIT. Axes-orthogonal splitting of worst LM into two partitions. The LM is splitted
% at its center coordinates.
%
%
% obj = LMSplit(obj, worstLM, splitDimension)
%
%
% INPUT
%
%   obj              object   HILOMOT object containing all relevant properties and methods.
%   worstLM:         1 x 1    Index of LM with the largest local loss function value.
%   splitDimension:  1 x 1    Index for the orthogonal splitting input dimension.
%
%
% 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

% Define sigmoid parameters for orthogonal split
if ~ischar(splitDimension)
    w = zeros(size(obj.zRegressor,2)+1,1);
    w(splitDimension) = 1;                                      % Split in dimension splitDimension
    w(end) = -obj.localModels(worstLM).center(splitDimension);  % Split at the center coordinate of the worst LM 
else
    w =  obj.localModels(worstLM).splittingParameter;
    w(end) = -obj.localModels(worstLM).center*obj.localModels(worstLM).splittingParameter(1:end-1);
end

% Perform split
[obj, forbiddenSplit] = LMPerformSplit(obj, w, worstLM);

% Check if split is forbidden, return without update
if forbiddenSplit
    return
end

% Evaluate new loss function values and update structure
idxLeafModels = find(obj.leafModels);
localLossFunctionValue = calcLocalLossFunction(obj, obj.unscaledOutput, obj.outputModel , obj.phi(:,idxLeafModels), obj.dataWeighting);
for k = 1:numel(idxLeafModels)
    obj.localModels(idxLeafModels(k)).localLossFunctionValue = localLossFunctionValue(k);
end

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

end


%% Function to perform split
function [obj, forbiddenSplit] = LMPerformSplit(obj, w, worstLM)

%
% 1) Split worst LM with sharp transitions
%

% Initialize flag forbidden split
forbiddenSplit = false;

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

% Get some constants
numberOfSamples = size(obj.input,1);
zLowerBound = min(obj.zRegressor);
zUpperBound = max(obj.zRegressor);

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

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

% Normalize and correct weight vector
deltaCenter = 0.01*(zUpperBound - zLowerBound);  % Choose a small value (sharp transition)
if any(deltaCenter<eps); deltaCenter=deltaCenter+eps; end
wNorm       = 20*w/(norm(w)*norm(deltaCenter));  % Normalize and correct weight vector

% Add two new children knot, i.e. local models, to the tree
obj = obj.addChildren(worstLM, wNorm);
obj = obj.addChildren(worstLM, -wNorm);

% Smoothness equal for all local models 
obj.localModels(end-1).localSmoothness = obj.smoothness;
obj.localModels(end).localSmoothness   = obj.smoothness;

% Calculate validity function values
obj.phi = [obj.phi zeros(numberOfSamples,2)];
psi     = 1./(1+exp((1./obj.smoothness) .* ( obj.zRegressor*wNorm(1:end-1) + wNorm(end)))); % Variante 2 (schneller)
psiComp = 1-psi;
obj.phi(:,end-1) = obj.phi(:,worstLM).*psi;
obj.phi(:,end)   = obj.phi(:,worstLM).*psiComp;

%
% 2) Calculate new LM centers
%

% Calculate centers from validity functions and data distribution (crisp transitions)
for dim = 1:size(obj.zRegressor,2)
   obj.localModels(end-1).center(dim) = obj.zRegressor(:,dim)'*obj.phi(:,end-1)/(sum(obj.phi(:,end-1))+eps);
   obj.localModels(end).center(dim)   = obj.zRegressor(:,dim)'*obj.phi(:,end)/(sum(obj.phi(:,end))+eps);
end

%
% 3) Correct adjustment of interpolation smoothness with updated centers
%

% Update phi with center info (correctly smoothed transition)
deltaCenter = obj.localModels(end-1).center - obj.localModels(end).center;    % Updated distance between LM-centers
if any(deltaCenter<eps); deltaCenter=deltaCenter+eps; end
wNorm       = 20*w/(norm(w)*norm(deltaCenter));                               % Normalize and correct weight vector
obj.localModels(end-1).splittingParameter = wNorm;                            % Update weights in LM objects
obj.localModels(end).splittingParameter   = -wNorm;

% Calculate phi values
psi     = 1./(1+exp((1./obj.smoothness) .* ( obj.zRegressor*wNorm(1:end-1) + wNorm(end)))); % Variante 2 (schneller)
psiComp = 1-psi;
obj.phi(:,end-1) = obj.phi(:,worstLM).*psi;
obj.phi(:,end)   = obj.phi(:,worstLM).*psiComp;

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


%
% 4) Estimate new LMs
%

% Estimate parameters of the LM
[obj.localModels(end-1).parameter, obj.localModels(end-1).pseudoInv] = estimateParametersLocal(obj,obj.xRegressor, obj.output, obj.phi(:,end-1), [], obj.dataWeighting);
[obj.localModels(end).parameter, obj.localModels(end).pseudoInv]     = estimateParametersLocal(obj,obj.xRegressor, obj.output, obj.phi(:,end), [], obj.dataWeighting);

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

end