function obj = LMSplitPostOpt(obj, objBest, worstLM)
% LMSPLITPOSTOPT. Update lastly added split of worst LM by applying a posterior split optimization.
%
%
% obj = LMSplitPostOpt(obj, objBest, worstLM)
%
%
% INPUT
%
%   obj       object   HILOMOT object containing all relevant properties and methods.
%   objBest:  object   Object containing the LMN with the best axes-orthogonal split as
%                      initialization.
%   worstLM:  1 x 1    Index of LM with the largest local loss function value.
%
%
% OUTPUT
%
%   obj       object   HILOMOT object containing all relevant properties and methods.
%

% HILOMOT Nonlinear System Identification Toolbox
% Benjamin Hartmann, 25-November-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2011 by Prof. Dr.-Ing. Oliver Nelles


%% Initialization

% Options for nonlinear optimization
% myOptions = optimset(...
%     'Display',            'off',...    % 'iter', 'final', 'off'
%     'Diagnostics',        'off',...
%     'MaxFunEvals',        1e3,...
%     'MaxIter',            1e3,...
%     'HessUpdate',         'bfgs',...   % 'bgfs', 'dfp', ('steepdesc')
%     'TolFun',             1e-10,...
%     'TolX',               1e-10,...
%     'FinDiffType',        'forward',...
%     'DiffMinChange',      1e-10,...
%     'LargeScale',         'off');

myOptions = optimset;
myOptions = optimset(myOptions,'Display','off');
% myOptions = optimset(myOptions,'Display','iter','Diagnostics','on');

% Switch of warnings
warning off optim:fminunc:SwitchingMethod
warning off MATLAB:rankDeficientMatrix

% Get initial splitting parameters from best orthogonal splitcl

wIni = objBest.localModels(end-1).splittingParameter;

% w0 is kept constant; w1, ..., wnz are going to be optimized
wRedIni = wIni(1:end-1);
w0 = wIni(end);

% Deliver weighted output of worst LM for better optimization performance
%weightedOutputWorstLM = obj.calcYhat(obj.xRegressor, obj.phi(:,worstLM), {obj.localModels(worstLM).parameter});

% Deliver output of the two lastly added LMs
weightedOutputWorstLM = obj.calcYhat(objBest.xRegressor, objBest.phi(:,end-1:end), {objBest.localModels(end-1:end).parameter});


% Index of worst LM from previous splitting
worstLM = objBest.history.splitLM(end);


%% Perform Optimization

if isempty(obj.scaleOutput)
    outputModelWithScaling = obj.outputModel;
else
    outputModelWithScaling = obj.scaleOutput.scale(obj.outputModel);
end

% Nonlinear split direction optimization
if ~obj.optGrad
    % Optimization with numerical approximation of gradient.
    %   Line of action:
    %     i)  First of all, try optimization with numerical gradient.
    %     ii) If i) fails, set parameters to initial values.
    try
        optGrad = 0;
        wRedOpt = fminunc(@(wRed) obj.obliqueGlobalLossFunction(wRed, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, weightedOutputWorstLM, ...
            obj.phi(:,worstLM), obj.smoothness, optGrad, obj.dataWeighting), wRedIni, myOptions);
    catch
        warning('hilomot:LMSplitOpt','Optimization aborted. Parameters are set to initial values!')
        wRedOpt = wRedIni;
    end
else
    % Optimization with analytical calculation of gradient.
    %   Line of action:
    %     i)   First of all, try optimization with analytical gradient.
    %     ii)  If i) fails, try optimization with numerical gradient.
    %     iii) If ii) fails, set parameters to initial values.
    try
        optGrad = 1;
        myOptions = optimset(myOptions,'GradObj','on','LargeScale','off');
        wRedOpt = fminunc(@(wRed) obj.obliqueGlobalLossFunction(wRed, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, weightedOutputWorstLM, ...
            obj.phi(:,worstLM), obj.smoothness, optGrad, obj.dataWeighting), wRedIni, myOptions);
    catch
        warning('hilomot:LMSplitOpt','Optimization with anlaytical gradient failed. Switch to numerical gradient.')
        try
            optGrad = 0;
            myOptions = optimset(myOptions,'GradObj','off');
            wRedOpt = fminunc(@(wRed) obj.obliqueGlobalLossFunction(wRed, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, weightedOutputWorstLM, ...
                obj.phi(:,worstLM), obj.smoothness, optGrad, obj.dataWeighting), wRedIni, myOptions);
        catch
            warning('hilomot:LMSplitOpt','Optimization aborted. Parameters are set to initial values (axes-orthogonal split)!')
            wRedOpt = wRedIni;
        end
    end
end
% Calculate loss function and model updates with optimal sigmoid parameters
[~, ~, parameterNew, centerNew, phiNew, outputModelNew] = obj.obliqueGlobalLossFunction(wRedOpt, w0, obj.xRegressor, obj.zRegressor, obj.output, ...
    outputModelWithScaling, weightedOutputWorstLM, obj.phi(:,worstLM), obj.smoothness, 0, obj.dataWeighting);


%% Update Local Model Network with new split information

% 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 sigmoid parameters

% Calculate normalization factor
deltaCenter = centerNew(1,:) - centerNew(2,:);
w           = [wRedOpt; w0];
kappa       = 20/(norm(w)*norm(deltaCenter));    % Factor that normalizes the sigmoid parameters

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

% Update LM informations
obj.localModels(end-1).localSmoothness = obj.smoothness; % Smoothness equal for all local models
obj.localModels(end).localSmoothness   = obj.smoothness;
obj.localModels(end-1).parameter       = parameterNew{1};
obj.localModels(end).parameter         = parameterNew{2};
obj.localModels(end-1).center          = centerNew(1,:);
obj.localModels(end).center            = centerNew(2,:);

% Store validity functions of the two new LMs
obj.phi = [obj.phi phiNew];

% Store current model output
if isempty(obj.scaleOutput)
    obj.outputModel = outputModelNew;
else
    obj.outputModel = obj.scaleOutput.unscale(outputModelNew);
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

% Here, the global loss function is evaluated separately, since the user defined loss criterion could differ
% from NRMSE that is used for nonlinear optimization
obj.history.globalLossFunction(end+1) = calcGlobalLossFunction(obj, obj.unscaledOutput, obj.outputModel, obj.dataWeighting);

% Check for sufficient amount of data samples
obj.idxForbiddenLM = [obj.idxForbiddenLM false(1,2)];
if sum(obj.phi(:,end-1)) < obj.numberOfPoints;
    obj.idxForbiddenLM(end-1) = true;
end
if sum(obj.phi(:,end)) < obj.numberOfPoints;
    obj.idxForbiddenLM(end) = true;
end

end


