function [CVScore nEff LOOE nEffAllLM] = calcLOOError(obj)
% CALCLOOERROR calculates the leave-one-out cross validation error (only
% the linear parameters are newly estimated) and determines the effective
% number of parameters by calculation of trace(S).
%
%
% [CVScore, nEff, LOOE] = calcPenaltyLossFunction(obj)
%
% OUTPUTS:
%
%   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)
%   LOOE:                   (N x q)     Leave-one-out CV error on each sample
%
% 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


if ~isempty(obj.scaleOutput)
    error('lossfunction:calcLOOError','calcLOOError is not useable with scaled outputs')
end
    
[N nx] = size(obj.xRegressor);
y      = obj.output;
yM     = obj.outputModel;
M      = sum(obj.leafModels);
q      = size(y,2);
X      = obj.xRegressor;

if isprop(obj,'phi')
    % hilomot models
    phi    = obj.phi;
else
    % lolimot models
    phi = obj.calculateVFV(obj.MSFValue);
    phi = cell2mat(phi);
end
lmIdx  = find(obj.leafModels);
SiiAll = zeros(N,M);

% Calculate and sum up diagonal elements Sii of the smoothing matrix for each local model. Use
% stored pseudo-inverse of WLS estimation, if available.
for k = 1:sum(obj.leafModels)
    if isprop(obj,'applySubsetSelection') && obj.applySubsetSelection
        X     = obj.xRegressor(:,obj.localModels(lmIdx(k)).parameter~=0);
        nx    = size(X,2);
    end
    r2        = phi(:,lmIdx(k)).^(1/2);
    XW        = X.*r2(:,ones(1,nx)).^3; % equivalent to Q*X*sqrt(Q)
    if any(strcmp(superclasses(obj),'sigmoidGlobalModel'))
        pseudoInv = obj.localModels(lmIdx(k)).pseudoInv;
        if isempty(pseudoInv) || (isprop(obj,'applySubsetSelection') && obj.applySubsetSelection)
            pseudoInv = pinv(X.*r2(:,ones(1,nx)));
        end
    else
        pseudoInv = pinv(X.*r2(:,ones(1,nx)));
    end
    for m = 1:nx
        SiiAll(:,k) = SiiAll(:,k) + XW(:,m).*pseudoInv(m,:)';
    end
    
    % Equivalent alternative 1:
    %
    %     for k2=1:N
    %         SiiAll(k2,lmIdx(k)) = phi(k2,lmIdx(k))*X(k2,:)*pseudoInv(:,k2)*sqrt(phi(k2,lmIdx(k)));
    %     end
    
    % Equivalent alternative 2:
    % (not implemented)
    %
    % Q = diag(obj.phi(:,lmIdx(k)));
    % S = S+Q*X*pinv(X'*Q*X)*X'*Q;
    % Sii(:,lmIdx(k)) = diag(S);
    
end
Sii     = sum(SiiAll,2);
LOOE    = (y-yM)./(ones(N,q)-Sii(:,ones(1,q)));
CVScore = sum(LOOE.^2)/N; % equivalent to MSE

% Number of effective parameters:
% Caution: Usually, for local model networks nEff is calculated as nEff =
% trace(S*S) [see e.g. Murray-Smith et al.]. However, this is
% computationally very expensive. Therefore, it will be calculated as nEff
% = trace(S). Consider: trace(S) >= trace(S*S) [Hastie et al.]. This leads
% to an overestimated amount of parameters which is more secure in terms of
% complexity considerations.
nEff = sum(Sii);
nEffAllLM = sum(SiiAll,1)';


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

