%% This function is used to generate a DPM using collapsed Dirichlet priors (CDP).
%% In this function a GMM is created for each phoneme using CDP. Classification is completed
%% by finding the phoneme whose GMM yields the highest posterior.

%% inputData contains all of the feature vectors (# samples x # features) and labels, outputData is the name/location
%% of the saved output, seed_val takes any values but if set to [] will use a random seed, truncLevel controls
%% the truncation level for the algorithm, and numIterations controls how many times the algorithm is run.


function mainFunc_cdp_avg(inputData, outputData, seed_val, truncLevel, numIterations)
    profile clear;
    profile on;
	% vdpgm is Kurihara's toolkit and should be added to this directory once downloaded
    addpath('vdpgm'); addpath('../common');
    
    eval(['load ' inputData]);

    % 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
    error_evl=zeros(numIterations,1);    
    error_val=zeros(numIterations,1);
    conf_matrix=cell(numIterations,1);
    conf_matrix_val=cell(numIterations,1);
    
    nummix_per_label = zeros(numIterations, max(labels_train));
    for j = 1:numIterations
        % Set model options
        opts = mkopts_cdp(truncLevel);     % Necessary to set to use CDP 
        opts.get_q_of_z = 1;	% Generate cluster labels (always set to 1)
        
        % 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

        % Build the DPM for each phoneme
        result = cell(max(labels_train),1);
        gaussians = cell(max(labels_train),1);
        pi_k = zeros(max(labels_train),1);
        tic;    % Start counting training time 
        for a = 1:max(labels_train)
			% Find features for a given label to train the model
			features = feat_vec_train(find(labels_train == a), :);
            % In CH-M there is an instance where a phoneme label does not exist
            % in the training data. In that case, we don't build a model for it                
            if isempty(features)                                                               
                continue;
            end
            pi_k(a) = sum(labels_train == a) / length(feat_vec_train);
            result{a} = vdpgm(features', opts);
            nummix_per_label(j,a) = sum(result{a}.hp_posterior.Nc >1);  %Find any mixture component that has 1 or more samples attributed to it.

            % Find all non-zero mixtures (i.e. at least one sample belongs to the mixture)
            mixtures = find(result{a}.hp_posterior.Nc >= 1);

            % Initialize the means and variances of each cluster
            mu = zeros(length(mixtures), size(feat_vec_train,2));
            sigma = zeros(size(feat_vec_train,2), size(feat_vec_train,2), length(mixtures));

            % This portion builds a GMM for each phoneme given the parameters found from the DPM
            for i = 1:length(mixtures)
                mu(i,:) =  result{a}.hp_posterior.m(:,mixtures(i))';
                tmp_sigma = result{a}.hp_posterior.B{mixtures(i)} / result{a}.hp_posterior.eta(mixtures(i));  
                tmp_sigma = triu(tmp_sigma) +triu(tmp_sigma)' - diag(diag(tmp_sigma));
                sigma(:,:,i) = tmp_sigma;
            end

            mix_weights = result{a}.hp_posterior.Nc(mixtures);
            gaussians{a} = gmdistribution(mu,sigma, mix_weights);
        end
        training_times(j) = toc;	%Record the training time for each phoneme model

		% Decode
        [output_cluster_assignment_evl] = choose_cluster(gaussians, pi_k, feat_vec_evltest);
        [acc, conf_matrix] = accuracy(output_cluster_assignment_evl,labels_evltest,unique(labels_train));

        [output_cluster_assignment_val] = choose_cluster(gaussians, pi_k, feat_vec_devtest);
        [acc_val, conf_matrix_val] = accuracy(output_cluster_assignment_val,labels_devtest,unique(labels_train));

        error_evl(j) = 1 - acc;
        error_val(j) = 1 - acc_val;
    end
	% Find average error rates
    avg_error_evl=mean(error_evl);
    avg_error_val=mean(error_val);
	% Find the average number of misture components for each phoneme model
    avg_nummix_per_label = mean(mean(nummix_per_label, 1));
    stats= profile('info');
    profile off
    save(outputData, 'avg_nummix_per_label', 'nummix_per_label', 'stats', 'training_times', 'avg_error_evl', 'avg_error_val', 'error_evl', 'error_val', 'conf_matrix', 'conf_matrix_val');

end
    
