function [query] = DoEQuery(obj)
% DOEQUERY places space-filling a query point in the worst LM region, considering already existing points (target input).
%
%
% [query] = DoEQuery(obj)
%
%
% INPUT
%
%   obj:    object   HILOMOTDOE object containing all relevant properties and methods.
%
%
% OUTPUT
%
%   query:  ? x p    New point(s) for measurement (query points).
%
%
% HiLoMoTDoE - Nonlinear System Identification Toolbox for Design of Experiments
% Benjamin Hartmann, 05-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


switch obj.query
    
    case 'worstLM' % Generate new point in worst performing LM.
        
        % Smoothness adjustment. Re-estimate LM parameters afterwards.
        obj.LMNObject.smoothness = obj.smoothnessQuery;
        if ~(obj.smoothnessQuery==obj.smoothnessSplit)
            obj.LMNObject.phi = obj.LMNObject.calculatePhi(obj.LMNObject.zRegressor, 1:length(obj.LMNObject.localModels));
            for k = 1:length(obj.LMNObject.localModels)
                % Estimate parameters of the LM
                obj.LMNObject.localModels(k).parameter = estimateParametersLocal(obj.LMNObject, obj.LMNObject.xRegressor, ...
                    obj.LMNObject.output, obj.LMNObject.phi(:,k), [], obj.LMNObject.dataWeighting);
            end
            obj.LMNObject.outputModel = obj.LMNObject.calculateModelOutput;
        end
        
        % Evaluate local loss functions of LMs with sharp transitions
        if obj.LMNObject.LOOCV
            nEffi = calcNEffi(obj.LMNObject); % effective number of parameters for all LM
            Ni    = sum(obj.LMNObject.phi(:,obj.LMNObject.leafModels),1)'; % Number of data points per LM
            J     = obj.LMNObject.calcLocalLossFunction(obj.LMNObject.output, obj.LMNObject.outputModel, obj.LMNObject.phi, [], [], 'SE'); % Sum of squared errors for each LM
            J     = sqrt(J(obj.LMNObject.leafModels)./(Ni-nEffi+eps)); % Standard deviation for each local model
        else
            J = obj.LMNObject.calcLocalLossFunction(obj.LMNObject.output, obj.LMNObject.outputModel, obj.LMNObject.phi, [], [], 'DRSE');
            J = J(obj.LMNObject.leafModels);
        end
        [~, worstLM] = max(J);
        leafModels   = find(obj.LMNObject.leafModels);
        worstLMIdx   = leafModels(worstLM);
        query        = obj.DoEGenerateNewData(worstLMIdx);
        
    case 'variance'
        
        % Maximum model output variance as query criterion.
        [covY] = obj.LMNObject.calcErrorbar(obj.candidatePoints);
        [~, maxIdx]  = max(covY);
        query = obj.candidatePoints(maxIdx,:);
        
        
        
        
end
end


function nEffi = calcNEffi(obj)

[N nx] = size(obj.xRegressor);
M      = sum(obj.leafModels);
X      = obj.xRegressor;
phi    = obj.phi;
lmIdx  = find(obj.leafModels);
SiiAll = zeros(N,M);
for k = 1:sum(obj.leafModels)
    
    pseudoInv = obj.localModels(lmIdx(k)).pseudoInv;
    r2        = phi(:,lmIdx(k)).^(3/2);
    XW        = X.*r2(:,ones(1,nx));
    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:
    %
    % Q = diag(obj.phi(:,lmIdx(k)));
    % S = S+Q*X*pinv(X'*Q*X)*X'*Q;
    % Sii(:,lmIdx(k)) = diag(S);
    
end

% 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.
nEffi = sum(SiiAll,1)';

end







