function [penaltyLossFunction, numberOfAllParameters, CVScore, nEff, AICc] = calcPenaltyLossFunction(obj)
% CALCPENALTYLOSSFUNCTION calculates the penalty loss function for the current LMN object. The model
% complexity is assessed with the Generalized Prediction Error (GPE), see [John Moody: "Prediction
% Risk and Architecture Selection for Neural Network", 1994] for details. By default, the number of
% model parameters is approximated with the number of all consequents regressors. If the flag LOOCV
% is used, the number of effective parameters is calculated with nEff = trace(S). Furthermore, if a
% nonlinear optimization of the sigmoid parameters is performed, the amount of parameters is roughly
% considered with nz*(M-1). This penalizes the higher model flexibility compared to only
% heuristically defined splitting.
%
% [penaltyLossFunction, numberOfAllParameters, CVScore, nEff] = calcPenaltyLossFunction(obj)
%
% OUTPUTS:
%
%   penaltyLossFunction:    (1 x 1)     penalty loss function value
%   numberOfAllParameters:  (1 x 1)     Approximated sum of linear and nonlinear
%                                       parameters in the LMN as complexity measure
%   CVScore:                (1 x 1)     Leave-one-out cross validation score (nonlinear
%                                       sigmoid parameter influences neglected)
%   nEff:                   (1 x 1)     Number of effective parameters. The degrees of
%                                       freedom for the weighted least squares estimation 
%                                       are considered with: nEff = trace(S)
%
% INPUTS:
%
%   obj                     (object)    global model object containing all relevant
%                                       properties and methods
%
%
% SYMBOLS AND ABBREVIATIONS
%
% LM:  Local model
%
% p:   Number of inputs (physical inputs)
% q:   Number of outputs
% N:   Number of data samples
% M:   Number of LMs
% nx:  Number of regressors (x)
% nz:  Number of regressors (z)
%
%
% LMNtool - Local Model Network Toolbox
% Benjamin Hartmann, 04-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


% Get constants
[N nz] = size(obj.zRegressor);
if isprop(obj,'leafModels')
    M  = sum(obj.leafModels);
else
    M = 1;
end

% Get current MSE value
if isprop(obj,'outputModel')
    errorMSE = calcGlobalLossFunction(obj, obj.unscaledOutput, obj.outputModel, obj.dataWeighting, obj.outputWeighting, 'MSE');
else
    errorMSE = calcGlobalLossFunction(obj, obj.unscaledOutput, obj.calculateModelOutput, obj.dataWeighting, obj.outputWeighting, 'MSE');
end

% Determine the overall number of LM parameters (all parameter not equal to 0)
numberOfAllParameters = sum(sum(cell2mat({obj.localModels(obj.leafModels).parameter})~=0));
if obj.LOOCV
    % Perform a leave-one-out cross validation (only the linear parameters are newly estimated) 
    % and determine the effective number of parameters by calculation of trace(S).
    [CVScore nEff]  = calcLOOError(obj);
else
    nEff = numberOfAllParameters;
    CVScore = [];
end

% Consider the contribution of the nonlinear sigmoid parameters heuristically
nEffAIC = nEff;
if any(strcmp(superclasses(obj),'sigmoidGlobalModel'))
    if isprop(obj,'oblique') && obj.oblique % case hilomot with oblique splitting
        % # all parameters = # LM parameters + # nonlinear splitting parameters ( nz*(M-1) )
        numberOfAllParameters = numberOfAllParameters + nz*(M-1);
        nEffAIC = nEff + nz*(M-1);
    end
end

% Calculate AIC value. If LOOCV is used, this leads to the Generalized Cross Validation (GCV) error
% [see e.g. C. LOADER: Local Regression and Likelihood (1999)]
AIC = N*log(errorMSE) + 2*nEffAIC;

% Use the corrected AIC (AICc) criterion. For large N AICc -> AIC. For small sample sizes (~N/nEff < 40) a bias-correction term is added to the standard
% AIC which leads to the AIC_corrected criterion [see  e.g. BURNHAM, ANDERSON: Model Selection and Multimodel Inference (2002)].
AICc = AIC + 2*nEffAIC*(nEffAIC+1)/(N-1-nEffAIC+eps);
if nEffAIC > N-1
    AICc = inf;
end

% Normalize to global loss function type
MSEaic = exp(AICc/N);
switch obj.lossFunctionGlobal
    
    case 'MSE'      % mean-squared-error
        penaltyLossFunction = MSEaic;
        
    case 'RMSE'     % root-mean-squared-error
        penaltyLossFunction = sqrt(MSEaic);
        
    case 'NMSE'     % normalized-mean-squared-error
        outputDifference2 = bsxfun(@minus,obj.unscaledOutput,mean(obj.unscaledOutput,1)).^2;
        den = sum(sum(outputDifference2));
        if den<eps; den=eps; end
        penaltyLossFunction = N*MSEaic/den;
        
    case 'NRMSE'    % normalized-root-mean-squared-error
        outputDifference2 = bsxfun(@minus,obj.unscaledOutput,mean(obj.unscaledOutput,1)).^2;
        den = sum(sum(outputDifference2));
        if den<eps; den=eps; end
        penaltyLossFunction = sqrt(N*MSEaic/den);
        
    otherwise
        error('lossfunction:calcPenaltyLossFunction','This type of lossfunction is not implemented so far. Choose "MSE", "RMSE", "NMSE" or "NRMSE"!')
        
end

end




