function [obj, forbiddenSplit] = LMSplitOpt(obj, wIni, worstLM)
% LMSPLITOPT. Axes-oblique splitting of worst LM into two partitions using nonlinear optimization.
%
%
% obj = LMSplitOpt(obj, wIni, worstLM)
%
%
% INPUT
%
%   obj       object   HILOMOT object containing all relevant properties and methods.
%   wIni:     vector   Vector of the best initial split.
%   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, 04-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 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');

% if obj.kStepPrediction > 1 % error
%     error('hilomot:LMSplitOpt','Split optimization is not supported for parallel simulation (dynamic models). Set obj.oblique=0!')
%     % parallel simulation is not implemeted in obliqueGlobalLossFunction.
%     % The calculation of the model output only works for staic models so
%     % far.
% end

% Initialize flag forbidden split
forbiddenSplit = false;

myOptions = optimset;

% myOptions = optimset(myOptions,'Display','off');
% myOptions = optimset(myOptions,'Display','iter');
% myOptions = optimset(myOptions,'Display','iter','Diagnostics','on');
myOptions = optimset(myOptions,'Display','off','Diagnostics','off','Algorithm','levenberg-marquardt');

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

% Get initial splitting parameters from best orthogonal split
% 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});

% Pre-calculate smoothing matrices of all local models except the worst LM before optimization for better performance
if obj.optLOOCV && sum(obj.leafModels)>1
    idxLM = setxor(find(obj.leafModels),worstLM);
    [N, nx] = size(obj.xRegressor);
    SiiAll = zeros(N,sum(obj.leafModels)-1);
    for k = 1:sum(obj.leafModels)-1
        r2        = obj.phi(:,idxLM(k)).^(1/2);
        XW        = obj.xRegressor.*r2(:,ones(1,nx)).^3; % equivalent to Q*X*sqrt(Q)
        pseudoInv = obj.localModels(idxLM(k)).pseudoInv;
        if isempty(pseudoInv)
            pseudoInv = pinv(X.*r2(:,ones(1,nx)));
        end
        
        for m = 1:nx
            SiiAll(:,k) = SiiAll(:,k) + XW(:,m).*pseudoInv(m,:)';
        end
    end
else
    SiiAll = [];
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, obj.LOOCV, obj.optLOOCV, SiiAll), wRedIni, myOptions);
%        wRedOpt = lsqnonlin(@(wRed) obj.obliqueGlobalLossFunction(wRed, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, weightedOutputWorstLM, ...
%            obj.phi(:,worstLM), obj.smoothness, optGrad, obj.dataWeighting, obj.LOOCV, obj.optLOOCV, SiiAll), wRedIni, [],[],myOptions);

        wRedOpt = lsqnonlin(@(wRed) obj.obliqueGlobalLossFunction(...
            wRed, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, ...
            weightedOutputWorstLM, obj.phi(:,worstLM), obj.smoothness, optGrad, ...
            obj.dataWeighting, obj.LOOCV, obj.optLOOCV, SiiAll, ...
            obj.kStepPrediction, obj.xInputDelay, obj.zInputDelay, obj.xOutputDelay, obj.zOutputDelay,...
            obj.leafModels, obj.localModels, worstLM), 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, obj.LOOCV, obj.optLOOCV, SiiAll,...
            obj.kStepPrediction, obj.xInputDelay, obj.zInputDelay, obj.xOutputDelay, obj.zOutputDelay,...
            obj.leafModels, obj.localModels, worstLM), 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, obj.LOOCV, obj.optLOOCV, SiiAll, ...
                obj.kStepPrediction, obj.xInputDelay, obj.zInputDelay, obj.xOutputDelay, obj.zOutputDelay,...
            obj.leafModels, obj.localModels, worstLM), 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, pseudoInv] = ...
    obj.obliqueGlobalLossFunction(...
    wRedOpt, w0, obj.xRegressor, obj.zRegressor, obj.output, outputModelWithScaling, ...
    weightedOutputWorstLM, obj.phi(:,worstLM), obj.smoothness, 0, ...
     obj.dataWeighting, obj.LOOCV, obj.optLOOCV, SiiAll,...
     obj.kStepPrediction, obj.xInputDelay, obj.zInputDelay, obj.xOutputDelay, obj.zOutputDelay,...
            obj.leafModels, obj.localModels, worstLM);

% Check for sufficient amount of data samples
if obj.pointsPerLMFactor * size(obj.xRegressor,2) - sum(phiNew(:,1)) > eps
    forbiddenSplit = true;
    return
end
if obj.pointsPerLMFactor * size(obj.xRegressor,2) - sum(phiNew(:,2)) > eps
    forbiddenSplit = true;
    return
end
        
        
%% Update Local Model Network with new split information

% Update active models
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];

% Normalize sigmoid parameters

% Calculate normalization factor
deltaCenter = centerNew(1,:) - centerNew(2,:);
if any(deltaCenter<eps); deltaCenter=deltaCenter+eps; end
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,:);

if obj.LOOCV
    obj.localModels(end-1).pseudoInv       = pseudoInv{1};
    obj.localModels(end).pseudoInv         = pseudoInv{2};
end

% 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);

end


