function J = calcLocalLossFunction(obj, output, outputModel, ...
    validityFunctionValue, dataWeighting, outputWeighting, criterion)
% CALCLOCALLOSSFUNCTION calculates the local loss function value for each
% column given in "validityFunctionValue". It is possible to calculate only
% the lossfunction values of certain local models, because every column is
% the validity of one local model. With the optional input
% "criterion" it is possible to change the error criterion. Therewith it is
% possible to calculate another error value without changing the criterion
% stored in the net object.
%
% J = calcLocalLossFunction(obj, output, outputModel, ...
%    validityFunctionValue, dataWeighting, outputWeighting, criterion)
%
% OUTPUTS:
%
%   J:                      (1 x 1)     Global loss function value
%
% INPUTS:
%
%   obj                     (object)    global model object containing all relevant
%                                       net and data set information and variables
%
%   output:                 (N x q)     Output matrix
%
%   outputModel:            (N x q)     Model output matrix
%
%   validityFunctionValue   (N x M)     validity function matrix
%
%   dataWeighting:          (N x 1)     data weighting matrix
%
%   outputWeighting:        (q x 1)     output weighting vector
%
%   criterion:              (string)    global error criterion. If not given,
%                                       the one specified in the object will be
%                                       used
%
%
% 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


if ~exist('criterion','var')
    % Read loss function criterion of the LMN object
    criterion = obj.lossFunctionLocal;
end

% Set dataWeighting to default values, if not defined already
[numberOfSamples numberOfOutputs] = size(output);
if ~exist('dataWeighting','var')
    dataWeighting = ones(numberOfSamples,1);
end
if isempty(dataWeighting)
    dataWeighting = ones(numberOfSamples,1);
end

% Set outputWeighting to default values, if not defined
if ~exist('outputWeighting','var')
    outputWeighting = ones(numberOfOutputs, 1);
end
if isempty(outputWeighting)
    outputWeighting = ones(numberOfOutputs, 1);
end

% Transform weighting vectors to weighting matrices in order to avoid expensive matrix
% multiplications
dataWeightingMatrix = dataWeighting(:,ones(1,numberOfOutputs));
outputWeightingMatrix = outputWeighting(:,ones(1,numberOfSamples))';

% Calculate squared loss
error2 = (output-outputModel).^2;

% Loss function evaluation
switch criterion
    
    case 'SE',  J = sum(validityFunctionValue'*(outputWeightingMatrix.*error2.*dataWeightingMatrix),2);
        
    case 'RSE', J = sqrt(sum(validityFunctionValue'*(outputWeightingMatrix.*error2.*dataWeightingMatrix),2));
        
    case 'DRSE'
        % DRSE: Density Root Square Error
        % Used for HILOMOT utilized as a Design of Experiments (DoE) algorithm.
        % The RSE is normalized with the sum of validity functions that represent
        % a measure of the number of local model data points.
        den = sum(validityFunctionValue,1)';
        if den<eps; den=eps; end
        J = sqrt(sum(validityFunctionValue'*(outputWeightingMatrix.*error2.*dataWeightingMatrix),2)./den);
            
    case 'MISCLASS',
        % Used for classification algorithm
        if size(output,2)==1 % Only one output (0/1)
            J = validityFunctionValue'*(dataWeighting.*(output~=(0.5*(1+sign(outputModel-0.5)))));
        else % Multiple outputs (0 0 1 ....)
            [~, index1]=max(output,[],2);      % Which output is correct?
            [~, index2]=max(outputModel,[],2); % Which output is currently most active?
            J = validityFunctionValue'*(dataWeighting.*(index1~=index2));
        end
        
    otherwise
        error('lossfunction:calcLocalLossFunction','This type of lossfunction is not allowed: Choose "SE", "RSE", "DRSE" or "MISCLASS"!')
        
end
