%% Routine for calculation of the number of effective parameters
function [nEff] = calcNEff(obj)

% Get the leafs of the local model object
leaves = find(obj.leafModels);

% Get the number of data points
numberOfData = size(obj.input,1);

% Check which type of local models are used
if any(strcmp(superclasses(obj),'sigmoidGlobalModel'))
    % for sigmoids
    validities = obj.phi(:,obj.leafModels);
elseif any(strcmp(superclasses(obj),'gaussianOrthoGlobalModel'))
    % for gaussians
    MSFValue = obj.MSFValue(obj.leafModels);
    validities = obj.calculateVFV(MSFValue);
end

% Initialize the smoothness matrix S
S = zeros(numberOfData,numberOfData);

% % "safety keyboard" if the numbers of elements dismatch
% if numel(validities) ~= numel(leaves)
%     keyboard
% end

% Calculate the smoothness matrix for each local model
for j=1:length(leaves) % for-loop over all LMs
    % Get the index of the current leaf
    currLeaf = leaves(j);
    
    % Calculate weight-matrix of the current leaf
    % Qj = diag(validitys(:,currLeaf));
    % Use a trick for better performance
    q = validities(:,j);
    Xj = obj.xRegressor(:, (obj.localModels(currLeaf).parameter~=0) );
    Qj = q(:,ones(size(Xj,2),1));
    
    % faster way of calculation
    if ~obj.optLOOCV
        B1 = Qj.*Xj; % replaced slow Q*X
        A = Xj'*B1;
        B2 = pinv(A);
        B3 = B1'; % replaced slow X'*Q
        Sj = B1*B2*B3;
    else
        pseudoInv = obj.localModels(currLeaf).pseudoInv;
        r         = validities(:,j).^(3/2);
        XW        = Xj.*r(:,ones(size(Xj,2),1));
        
        Sj = XW*pseudoInv;
    end
    
    % Slow calculation method
    % Sj = (Qj*Xj)*pinv(Xj'*(Qj*Xj))*(Xj'*Qj); % old, slow calculation method
    
    S = S + Sj; % super-positioned global smoothing matrix
end

% Calculate the trace of the 
% nEff = trace(S'*S);
% Use a trick for better performance
nEff = sum(sum(S.^2));
 
end
