function [means, variances, weights, minLogL, Si, expectation, posteriors] = fitGMMwithPrior(...,
    data, ...
    prior, ...
    priorMap, ...
    means, ...
    variances,...
    weights, ...
    maxIters, ...
    convCriterion, ...
    numBins, ...
    verbosity, ...
    gmmFig, ...
    costFig)
% Description: Fits a Gaussian mixture model to data by means of Expectation Maximization
%
% Call: [means, variances, weights, minLogL, Si, dataEstimate] = fitGMM(data, means, variances, weights, useFixedMeans,
%           useSameVariance, maxMMIters, convCriterion, numBins, verbosity, figToPlotTo)
%
% Input:    data:               Data to fit the GMM to
%           prior:              Prior on tissue or structure classes
%           means:              Means of Gaussians in mixture
%           variances:          Variances of Gaussians in mixture
%           weights:            Weights of Gaussians in mixture
%           usedFixedMeans:     Means are updated according to the min and max value of the data (hence 2 parameters)
%           useSameVariance:    Equal variance is computed for all Gaussians
%           maxIters:           Maximum number of iterations allowed for fit
%           convCriterion:      Convergence threshold (relative change in cost)
%           numBins:            bins to plot histogram with
%           verbosity:          Verbosity level, higher value provides more output to stdout
%           figToPlotTo:        Figure handle to plot mixture model fitting to as iterations progress
%
% Output:   means:              Gaussian means that optimize likelihood of data
%           variances:          Gaussian variances that optimize likelihood of data
%           weights:            Gaussian weights that optimize likelihood of data
%           minLogL:            The negative loglikelihood at convergence
%           Si:                 The weights (or confidence) of each voxels in fit (posterior probability over variance)
%           dataEstimate:       The estimated data, given the optimized parameters
%
% Author: Christian Thode Larsen, christian.thode.larsen@gmail.com

% optional parameter initialization
if nargin < 7
    maxIters = 100;
end
if nargin < 8
    convCriterion = 1e-5;
end
if nargin < 9
    numBins = 1024;
end
if nargin < 10
    verbosity = 'none';
end
if nargin < 11
    gmmFig = -1;
end
if nargin < 12
    costFig = -1;
end

% initialize
numGaussians = length(means);
numClasses = size(prior, 2);
numVoxels = numel(data);
squaredDists = zeros(numVoxels, numGaussians);
likelihoods = zeros(numVoxels, numGaussians);

% history of minLogLikelihoods - used to determine convergenge between iterations
costHist = 1/eps;
converged = false;

if ishandle(gmmFig)
    gaussianEvals = zeros(numBins, numGaussians);
    classEvals = zeros(numBins, numClasses);
    
    [histogram, binCenters] = hist(data, numBins);
    histogram = histogram / sum(histogram);
    binWidth = binCenters(2) - binCenters(1);

    priorFrequency = sum(prior) / numVoxels;
end

% Performing EM
for mmIter = 1:maxIters
    % plot the state of affairs
    if ishandle(gmmFig)
        figure(gmmFig)
        for i = 1:numGaussians
            gaussianEvals(:, i) =  (1 / sqrt(2*pi*variances(i))) * ...
                exp(-0.5 * (binCenters - means(i)).^2 / variances(i)) * weights(i) * priorFrequency(priorMap(i));
        end
        gaussianEvals = gaussianEvals * binWidth;
        for i = 1:numClasses
            classEvals(:,i) = sum(gaussianEvals(:,i == priorMap),2);
        end
        clf
        grid
        bar(binCenters, histogram);
        hold on
        plot(binCenters, classEvals', 'lineWidth', 1);
        plot(binCenters, sum(gaussianEvals, 2), 'r', 'lineWidth', 1);
        hold off
        drawnow
    end
    
    for i = 1:numGaussians
        squaredDists(:, i) = (data - means(i)).^2;
        likelihoods(:, i) = (1 / sqrt(2*pi*variances(i))) * exp(-0.5 * squaredDists(:,i) / variances(i)) * weights(i);
    end
    likelihoods = likelihoods .* prior(:,priorMap);
    
    normalizer = sum(likelihoods, 2);
    
    minLogL = -sum(log(normalizer + eps));
    costHist(end + 1) = minLogL;%#ok
    
    relChange = (costHist(end - 1) - costHist(end)) / numVoxels;
    if abs(relChange) < convCriterion
        if ~isempty(strfind(verbosity, 'debug'))
            fprintf('\tGGM fit converged in %i iterations\n', mmIter);
        end
        converged = 1;
        break
    end
    
    if ishandle(costFig)
        figure(costFig)
        plot(costHist(2:end));
        axis tight
        title('Cost');
        xlabel('Iteration');
        ylabel('Cost');        
        drawnow
    end
    
    posteriors = likelihoods ./ repmat(normalizer + eps, [1, numGaussians]);
    posteriorsSum = sum(posteriors);
    posteriorsSumSum = zeros(numClasses, 1);
    
    for g = 1:numGaussians  
        l = priorMap(g);
        posteriorsSumSum(l) = posteriorsSumSum(l) + sum(posteriorsSum(g));
    end

    % update weights
    for g = 1:numGaussians
        l = priorMap(g);
        weights(g) = posteriorsSum(g) / posteriorsSumSum(l);
    end

    % verify that \sum w_gl = 1
    for l = 1:numClasses
        if(abs(sum(weights(priorMap == l)) - 1) > 1e-7)
            error('Weights do not sum up to one')
        end
    end
    
    % update means
    means = sum(posteriors .* repmat(data,[1, numGaussians])) ./ (posteriorsSum + eps);
    
    % update variances
    for i = 1:numGaussians
        squaredDists(:, i) = (data - means(i)).^2;
    end 
    variances = sum(posteriors .* squaredDists) ./ (posteriorsSum + eps);
    
end

if (~converged && ~isempty(strfind(verbosity, 'debug')))
    disp(['GMM stopped: Maximum number (', num2str(maxIters), ') of iterations reached'])
end

if(nargout > 4)
    Sik = posteriors ./ repmat(variances, [numVoxels, 1]);
    Si = sum(Sik, 2);
end

if(nargout > 5)
    expectation = sum(Sik .* repmat(means, [numVoxels, 1]), 2) ./ (Si + eps);
end



end
