function [AllLMN alpha YMtrain trainRMSE] = modelByCommittee(data,testData)

% data:      Training data as [u1,u2,...,up y]

% Select training algorithms
% methods = {'hilomot','hilomotQuad','hilomotCube','hilomot+','hilomot+Quad','hilomot+Cube'};
% methods = {'hilomot','hilomotQuad','hilomotCube','hilomot+','hilomot+ x2','hilomot+ x4'};
% methods = {'hilomotCube','hilomot+ x4'};
methods = {'hilomotQuad','hilomotCube','hilomot+Quad','hilomot+Cube'};
% methods = {'hilomot+Quad'};

% Initialization
nCom = length(methods);
N    = size(data,1);
p    = size(data,2)-1; 
AllLMN   = cell(2,nCom);
YM       = zeros(N,nCom);
LOOError = zeros(1,nCom);



for idxMethod = 1:nCom
    
    switch methods{idxMethod}
        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;
        case 'hilomotQuad'
            LMN = hilomot;
            LMN.xRegressorDegree = 2;
        case 'hilomotCube'
            LMN = hilomot;
            LMN.xRegressorDegree = 3;
        case 'hilomotSparseQuad'
            LMN = hilomot;
            LMN.xRegressorDegree = 2;
            LMN.xRegressorType = 'sparsePolynomial';
        case 'hilomot+'
            LMN = hilomotSubset;
        case 'hilomot+ x0.5'
            LMN = hilomotSubset;
            LMN.increasingNumberOfLMParameters = size(data,2)/2;
        case 'hilomot+ x2'
            LMN = hilomotSubset;
            LMN.increasingNumberOfLMParameters = 2*size(data,2);
        case 'hilomot+ x4'
            LMN = hilomotSubset;
            LMN.increasingNumberOfLMParameters = 4*size(data,2);
        case 'hilomot+Quad'
            LMN = hilomotSubset;
            exponentMatrix = LMN.buildExponentMatrix('polynomial',2,p);
            LMN.increasingNumberOfLMParameters = size(exponentMatrix,1);
        case 'hilomot+Cube'
            LMN = hilomotSubset;
            exponentMatrix = LMN.buildExponentMatrix('polynomial',3,p);
            LMN.increasingNumberOfLMParameters = size(exponentMatrix,1);
        otherwise
            error('LMNTrain:unknownMethod','Method "%s" is unknown.',methods{idxMethod})
    end
    if isprop(LMN,'optLOOCV')
        LMN.optLOOCV = false;
    end
    if strncmp(methods{idxMethod},'hilomot+',8)
        LMN.xRegressorDegree = 3;
        LMN.estimationProcedure = 'REDUCED_STEPWISE_LS';
    end



    
    
    % Maximal training time per algorithm in minutes
    LMN.history.displayMode = true;
    LMN.lossFunctionGlobal = 'RMSE';
    LMN.smoothness = 1;
    LMN.LOOCV = true;
    LMN.maxTrainTime = inf;
    % 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);
    if nargin>1
        LMN.testInput  = testData(:,1:end-1);
        LMN.testOutput = testData(:,end);
    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.outputModel;
        
        % 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})
%         keyboard
%     end
    
    % Clearing LMN
    clear LMN
    
end

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

% Evaluate model on training data
if nargout > 2
    YMtrain = YM*alpha';
end
% Calculate error on training data
if nargout > 3
    trainRMSE = AllLMN{2,1}.calcGlobalLossFunction(data(:,end),YMtrain);
end

