function [varargout] = aamlib_gaussianClassifier(varargin)
% function [varargout] = ammlib_groupClassifier(varargin)
%
% Dr. A. I. Hanna (2007)

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename,' must come param/value in pairs.'])
end
groups = [];
pc_ind = [];
X = [];
% just to make sure that the covariance matrix is not singular
reg_const = 0.1;
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'groups'
            groups = varargin{i+1};
        case 'data'
            X = varargin{i+1};
        case 'reg_const'
            reg_const = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(groups)
    error('Parameter groups must be supplied');
end
if isempty(X)
    error('Parameter X must be supplied');
end
N = size(X,2);
G = length(groups);
gaussProbDstns = cell(G,1);
for i=1:G
    b = [groups(i).shapedata(:).b];
    gaussDist.mean = mean(b,2);
    bc = b - gaussDist.mean*ones(1, size(b,2));
    gaussDist.sigma = cov(bc')+ eye(size(bc,1))*reg_const;
    gaussProbDstns{i} = gaussDist;
end
g = [];
for j=1:G
    mu = gaussProbDstns{j}.mean;
    g = cat(1, g, multiGauss(X, mu, gaussProbDstns{j}.sigma));
end
[g, ind] = sort(g, 1, 'descend');
varargout{1} = ind;
varargout{2} = g;
return;
%%%
%
%%%
function g = multiGauss(X, mu, sigma)
N = length(X);
inv_sigma = pinv(sigma);
sqrt_det_sigma = sqrt(det(sigma));
c = (2*pi)^(N/2)* (sqrt_det_sigma + realmin);
g = zeros(1, size(X,2));
for i = 1:size(X,2)
    x = X(:,i);
    x = x - mu;
    e = exp(-0.5*(x'*inv_sigma*x));
    g(i) = e./c;
end
return