function J = calcGlobalLossFunction(obj, output, outputModel, dataWeighting, outputWeighting, criterion)
% CALCGLOBALLOSSFUNCTION calculates the global loss function for the net
% object. 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 which is stored in the net object.
%
% J = calcGlobalLossFunction(obj, output, outputModel, 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
%
%   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
% Tobias Ebert, 17-Nov-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles

% 18.11.11  help updated (TE)
% 16.12.11  Added "eps" to outputDifference2, otherwise J=NaN possible (B.Hartmann)


if ~exist('criterion','var')
    % Used criterion of the lolimot object
    criterion = obj.lossFunctionGlobal;
end

% Get some constants
[numberOfSamples numberOfOutputs] = size(output);


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

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

% assemble weightings to matrix to avoid matrix multiplication
dataWeightingMatrix = dataWeighting(:,ones(1,numberOfOutputs));
outputWeightingMatrix = outputWeighting(:,ones(1,numberOfSamples))';

error2 = (output-outputModel).^2;   % N x q squared error matrix

% Loss function evaluation
switch criterion
    
    case 'MSE'      % mean-squared-error
        J = mean(sum(outputWeightingMatrix.*error2.*dataWeightingMatrix,2));
        
    case 'RMSE'     % root-mean-squared-error
        J = sqrt(mean(sum(outputWeightingMatrix.*error2.*dataWeightingMatrix,2)));
        
    case 'NMSE'     % normalized-mean-squared-error
        outputMean = mean(output,1);
        outputDifference2 = (output-outputMean(ones(numberOfSamples,1),:)).^2;
        den = sum(sum(outputDifference2,2));
        if den<eps; den=eps; end
        J = sum(sum(outputWeightingMatrix.*error2.*dataWeightingMatrix,2))/den;
        
    case 'NRMSE'    % normalized-root-mean-squared-error
        outputMean = mean(output,1);
        outputDifference2 = (output-outputMean(ones(numberOfSamples,1),:)).^2;
        den = sum(sum(outputDifference2,2));
        if den<eps; den=eps; end
        J = sqrt(sum(sum(outputWeightingMatrix.*error2.*dataWeightingMatrix,2))/den);
        
    case 'MISCLASS' % Used for classification algorithm
        if size(output,2)==1 % Only one output (0/1)
            J = 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 = dataWeighting'*(index1~=index2);  % How many misclassifications?
        end
        
    case 'R2'       % coefficient of determination aka R squared
%         outputMean = mean(output,1);
%         varY = sum((output-outputMean).^2,1);
%         numerator = sum((outputModel-outputMean).^2,1);
%         if varY<eps; varY=eps; end
%         J = numerator/varY;        
        outputMean = mean(output,1);
        varY       = sum((output-outputMean).^2);
        
%         SSE=sum((y-ypred).^2);
%         SST=sum((y-mean(y)).^2);
%         
%         SSE
%         SSE =
%         0.166452761430755
%         
%         SST
%         SST =
%         1.10029991380876
%         
%         1- SSE/SST
%         ans =
        
        
        if varY<eps; varY=eps; end
        J = 1 - sum(error2)/varY;
        
    otherwise
        error('lossfunction:calcGlobalLossFunction','This type of lossfunction is not allowed. Choose "MSE", "RMSE", "NMSE", "NRMSE", "MISCLASS" or "R2"!')
        
end
