% This function is used for TIMIT data that was generated by an HTK alignment (not the manual transcriptions)

% This function generates GMMs with k mixtures from the trainining data 
% (inputData), and saves it to outputData. 

%% inputData contains all of the feature vectors (# samples x # features) and labels, k is the number
%% of mixture components to use for each phoneme label, outputData is the name/location
%% of the saved output,  numIterations controls how many times the algorithm is run (in this case the best
%% of these iterations is returned, not the average), and seed_val can be used to control the seeding of the
%% initialization.

function mainFunc_htk(inputData, k, outputData, numIterations, seed_val)
    addpath('../common');	% Incude basic functions in the path
    eval(['load ' inputData]);	%load the data
    
    % This sets the state of the random initialization for the model
    if ~isempty(seed_val)
        opts.seed = seed_val;
    else
        system('sleep 1');
        rng('default'); rng('shuffle');     %You have to reinitialize the state of rand to use 'shuffle' (clock based)
        opts.seed = randi(2^32);
    end

    [n_examples, n_features] = size(feat_vec_train);

    % Normalize the data
    [meanv, stdv, feat_vec_train] = normalize(feat_vec_train, [], []);

    % Use the means and variance from the training to normalize the validation 
    % and test data
    [~, ~, feat_vec_devtest] = normalize(feat_vec_devtest, meanv, stdv); 
    [~, ~, feat_vec_evltest] = normalize(feat_vec_evltest, meanv, stdv); 

    % Initialize the GMMs
    gaussians = cell(max(labels_train), 1);
    pi_k = zeros(max(labels_train),1);

    % Set options for the GMM training
    options = statset('Display','final', 'MaxIter', 300);
    bad_labels = zeros(max(labels_train),1);

    tic; %Measure training time
    for i = 1:max(labels_train)
        if sum(labels_train == i) < n_features
            bad_labels(i) = 1;
            continue;
        end
		
        % Find all instances of a given phoneme
        features = feat_vec_train(find(labels_train == i), :);
        pi_k(i) = sum(labels_train == i) / length(feat_vec_train);
        
        % We assume the covariance is diagonal, we add a small number to each
        % index of the covariance matrix to avoid ill-conditioned matrix
        % errors, and the 'Replicates' option randomly initializes the model
        % numIterations times to find the best Gaussian model

        % Train a GMM for each phoneme
        gaussians{i} = gmdistribution.fit(features,k,'Options',options, 'CovType', 'diagonal', 'Regularize', 1e-15, 'Start', 'randSample', 'Replicates', numIterations); 
        %gaussians{i} = gmdistribution.fit(features,k,'Options',options, 'Regularize', 1e-15, 'Start', 'randSample', 'Replicates', numIterations); 
    end
    training_time = toc;	%Record training times
    avg_training_time=training_time/numIterations;

	% Decode
    [output_evl_labels] = choose_cluster(gaussians, pi_k, feat_vec_evltest);
    [output_val_labels] = choose_cluster(gaussians, pi_k, feat_vec_devtest);
    
    % Determine misclassification error on the validation set
    [acc_val, conf_matrix_val] = accuracy(output_val_labels,labels_devtest,unique(labels_train));
   
    % Determine misclassification error on the evaluation set
    [acc_evl, conf_matrix_evl] = accuracy(output_evl_labels,labels_evltest,unique(labels_train));

    clear feat_vec_train feat_vec_devtest feat_vec_evltest labels_train labels_devtest labels_evltest;
    save(outputData);
end
