function [G,Z] = gmmemdistances(D)
% GMMEMDISTANCES Fit a Gaussian mixture model to a sample of distances
%
%  [G,Z] = GMMEMDISTANCES(D) where D is a vector containing a
%  sample of distances. The result G is the Netlab GMM struct and Z
%  are the posteriors of the data under the model G.

% Copyright (c) 2012
% Jerod Weinman
% jerod@acm.org

% $Id: gmmemdistances.m 595 2014-10-24 18:15:54Z weinman $
 
% Jerod Weinman, Zach Butler, Dugan Knoll, Alexander Turner, Yirong Jing
% jerod@acm.org
% Copyright 2011, 2012, 2013, 2014
 
%    This file is part of PIONEER WordNorm.
%
%    PIONEER WordNorm is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    PIONEER WordNorm is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with PIONEER WordNorm.  If not, see <http://www.gnu.org/licenses/>.

% Number of components
K=min(4,length(unique(D)));

% Initialize by k-means
opt = zeros(1,14);
opt(1) = -1; % turn display off completely
%  opt(2) = kmeans_tol; % For K-Means
opt(3) = 1e-8; % For GMM

opt(14) = 100; % iter limit

Ds = sort(D);

while K>0
    
    % Initial mixture model
    G0 = gmm(1,K,'diag');

    % Initialize centres to percentiles, trimming ends for stability to outliers 
    P = linspace(0,1,K+2); 
    G0.centres = Ds(ceil(P(2:end-1)*length(Ds)));
    
    
    % Initialize EM structure
    opt(5) = 0; % do not randomly initialize centers
    G = gmminit(G0, D, opt);
    
    % Run EM 
    opt(5) = 1; % check for min covars
    [G,opt] = gmmem(G, D, opt);
    
    if ~any(isnan(G.centres))
      break
    else
      K = K-1;
    end;
end;

if K>0
  % Get component posteriors
  Z = gmmpost(G,D);
else
  Z = ones(size(D));
  G = gmm(1,1,'diag');
  G.centres = mean(D);
  G.var = var(D);
end
