% This function is used to find the error generated by random forest models for
% data that is already partitioned into training, validation, and test
% sets. inputData is the source data containing the feature vectors where
% the # of rows represent the # of samples and columns represent # of features
% (e.g. 'feat_vec_foo' - foo = 'train', 'val', or 'evl'), and labels (e.g. 
% 'labels_foo'), numTrees is the number of trees to use in the ensemble,
% outputData is the name of the .mat file to save, and max iter determines
% how many iterations to use for averaging error rates

function mainFunc(inputData, numTrees, outputData, maxiter)
    % Add basic functions to path
    addpath('../common');
    eval(['load ' inputData]);

    [n_examples, n_features] = size(feat_vec_train);

    error = zeros(1,maxiter);
    error_val = zeros(1,maxiter);
    confm=zeros(length(unique(labels_train)),length(unique(labels_train)));   
    disp(['Started at: ' mat2str(clock)]);
    tr = [feat_vec_train labels_train];
    val = [feat_vec_devtest labels_devtest];
    test = [feat_vec_evltest labels_evltest];

    % Randomly mix the data up
    p1 = randperm(size(tr,1));
    p2 = randperm(size(val,1));
    p3 = randperm(size(test,1));
    tr = tr(p1,:);
    val = val(p2,:);
    test = test(p3,:);
    
    % Normalize all the features to mean zero and variance 1 using the 
    % training data.
    [meanv, stdv, tr(:, 1 : n_features)] = normalize(tr(:, 1 : n_features), [], []);
    [~, ~, val(:, 1 : n_features)] = normalize(val(:, 1 : n_features), meanv, stdv); 
    [~, ~, test(:, 1 : n_features)] = normalize(test(:, 1 : n_features), meanv, stdv); 

    training_times=zeros(maxiter,1);
    for itr = 1:maxiter
        % Train a random forest and then use compact to convert it into a
        % manageable size to save
    	tic;    % Begin recording training time
        b = compact(TreeBagger(numTrees,tr(:, 1 : n_features),tr(:,n_features+1),'method','classification')); %, 'OOBPred','on'));
        training_times(itr)=toc;    % Measure training time
        prediction_val=cellfun(@str2num, predict(b,val(:, 1 : n_features)));
        prediction=cellfun(@str2num, predict(b,test(:, 1 : n_features)));
        
        % Calculate the accuracy
        [acc, conf_matrix] = accuracy(prediction,test(:,n_features+1),unique(labels_train));
        confm=confm+conf_matrix;
        error(itr)=1-acc;

        [acc_val, conf_matrix_val] = accuracy(prediction_val,val(:,n_features+1),unique(labels_train));
        error_val(itr)=1-acc_val;

        disp(['Iteration ' num2str(itr) ' finished at ' mat2str(clock)]);
        clear b;
    end
    
    % Find averages
    avg_training_time=sum(training_times)/maxiter;
    avg_error = sum(error)/maxiter;
    avg_confm = confm/sum(confm(:));
    avg_error_val = sum(error_val)/maxiter;

    disp(['Finished at: ' mat2str(clock)]);
    clear labels_train labels_devtest labels_evltest feat_vec_train ...
        feat_vec_evltest feat_vec_devtest tr val test p1 p2 p3;
    save(outputData, '-v7.3');
end
