function y = evalModelByCommittee(input, AllLMN, alpha)

% Calculate model output of a model committee given by AllLMN and alpha.


% data:      Training data as [u1,u2,...,up y]
% testInput: Inputs [u1,u2,...,up] to calculate model output YMtest


% Select training algorithms
methods = {'hilomot','hilomotQuad','hilomot+'};

% Initialization
AllLMN   = cell(1,length(methods));
YM       = zeros(size(testInput,1),length(methods));
LOOError = zeros(1,length(methods));

% Number of committee members
nCom = length(methods);

for idxMethod = 1:nCom
    
    switch methods{idxMethod}
        case 'globalPolynom'
            LMN = globalPolynom;
        case 'lolimot'
            LMN = lolimot;
        case 'lolimotQuad'
            LMN = lolimot;
            LMN.xRegressorDegree = 2;
        case 'lolimotSparseQuad'
            LMN = lolimot;
            LMN.xRegressorDegree = 2;
            LMN.xRegressorType = 'sparsePolynomial';
        case 'hilomot'
            LMN = hilomot;
            LMN.optLOOCV = true;
        case 'hilomotQuad'
            LMN = hilomot;
            LMN.xRegressorDegree = 2;
            LMN.optLOOCV = true;
        case 'hilomotSparseQuad'
            LMN = hilomot;
            LMN.xRegressorDegree = 2;
            LMN.xRegressorType = 'sparsePolynomial';
            LMN.optLOOCV = true;
        case 'lolimotMerge'
            LMN = lolimotMerge;
        case 'polymot'
            LMN = polymot;
            LMN.xRegressorDegree = 3;
        case 'polymotMerge'
            LMN = polymotMerge;
            LMN.xRegressorDegree = 3;
        case 'hilomot+'
            LMN = hilomotSubset;
            LMN.xRegressorDegree = 3;
        case 'hilomot2+'
            LMN = hilomotSubset;
            LMN.xRegressorDegree = 3;
            LMN.increasingNumberOfLMParameters = size(data,2)/2;
        case 'hilomot3+'
            LMN = hilomotSubset;
            LMN.xRegressorDegree = 3;
            LMN.increasingNumberOfLMParameters = 2*size(data,2);
        otherwise
            error('LMNTrain:unknownMethod','Method "%s" is unknown.',methods{idxMethod})
    end


    
    
    % Maximal training time per algorithm in minutes
    LMN.history.displayMode = false;
    LMN.lossFunctionGlobal = 'RMSE';
    LMN.smoothness = 1.25;
    LMN.LOOCV = true;
    LMN.maxTrainTime = 1;
    % LMN.maxNumberOfLM = 1;
    % LMN.scaleInput = [];
    % LMN.maxNumberOfParameters = 20;
    
    % Generate training data and store in hilomot object
    LMN.input  = data(:,1:end-1);
    LMN.output = data(:,end);
    
    
    
    % Display
    if LMN.history.displayMode
        fprintf(['\nCurrent training method: ',methods{idxMethod},'\n'])
    end
    
    % Perform training
    try
        LMN = LMN.train;
        LMN.leafModels = LMN.history.leafModelIter{LMN.suggestedNet};
        
        % Storage LMN and methods name in variable
        AllLMN{1,idxMethod} = methods{idxMethod};
        AllLMN{2,idxMethod} = LMN;
        
        YM(:,idxMethod) = LMN.calculateModelOutput(testInput);
        
        % Calculate LOO-error. MSE is equivalent to estimated model variance
        errLMN = LMN;
        errLMN.lossFunctionGlobal = 'MSE'; 
        LOOError(idxMethod) = errLMN.calcLOOError;
    catch
        warning('Training aborted!')
        disp(methods{idxMethod})
    end
    
    % Clearing LMN
    clear LMN
    
end

% Calculate weighting of committee members as normalized inverse model variance.
LOOErrorInv = 1./(LOOError+eps);
alpha       = LOOErrorInv/sum(LOOErrorInv+eps)

% Calculate variance of committee
YMmean = mean(YM,2);
YMmean = YMmean(:,ones(length(methods),1));
varVec = 1/(length(methods)-1)*sum((YM-YMmean).^2,2);

% Find point of biggest committee variance
[~, maxIdx]  = max(varVec);
query = testInput(maxIdx,:);

% Evaluate model for new data
N       = size(data,1);
NG      = size(testInput,1);
YMtrain = zeros(N,1);
YMtest  = zeros(NG,1);
for k = 1:nCom
    YMtrain = YMtrain + alpha(k)*AllLMN{2,k}.outputModel;
    YMtest  = YMtest  + alpha(k)*AllLMN{2,k}.calculateModelOutput(testInput);
end


trainRMSE = AllLMN{2,1}.calcGlobalLossFunction(data(:,end),YMtrain);


% if size(data,2)==2
%     figure(98)
%     clf
%     drawnow
%     hold all;
%         plot(data(:,1),data(:,2),'bx','markersize',12,'linewidth',2)
%     plot(u,y,'b-','linewidth',2)
% 
%     for i = 1:length(methods)
%         plot(candidatePoints,YM(:,i))
% %             pause
%     end
%     % plot(candidatePoints,mean(YM,2),'b--','linewidth',2)
%     
%     
% %     plot(candidatePoints,YMw+5*varVec,'r--','linewidth',1)
% %     plot(candidatePoints,YMw-5*varVec,'r--','linewidth',1)
%     %
%     plot(query,YMw(maxIdx),'ko','markersize',14,'linewidth',2)
%     
%     
%     plot(candidatePoints,YMw,'r','linewidth',2)
%     
% end
% pause%(0.5)

