function parm = init_gmix(data, n_modes, names, min_std, rand_init);
% function parm = init_gmix(data,n_modes,names,min_std,[rand_init]);
%
% INPUT:
% Stand-alone version of gmix_init (does not require toolkit)
% data: numFeature rows by numSamples columns. [n_feat x n_samp]
% n_modes: number of gaussian modes
% names: a cellarray of strings
% min_std: minimum standard deviation
% rand_init: Random initialization is on by default, but can be disabled by setting
%   rand_init equal to zero.  In that case, the first n_modes samples are
%   used to initialize the means.
%
% OUTPUT:
%   parm: TBD
%

%
% NAVAL UNDERSEA WARFARE CENTER
% class-specific@npt.nuwc.navy.mil
% http://www.npt.nuwc.navy.mil/csf
% Paul M. Baggenstoss
% 401-832-8240
% Do not disseminate.
%



%  Date              Reason
%  ------------      --------------
%  Mar 28, 1997      Initial release
%  May 21, 1997      Handles non-integer nmode by rounding
%  Dec 10, 1999      Upgraded to MATLAB 5
%  Dec 20, 1999      took out normalization of min_std
%  Feb 27, 2001      Changed input arguments to be compatible with toolkit.
%  Apr 3, 2001       Disabled data_means, data_std and added check for ill-conditioning.
%  Apr 9, 2004       CAS: fixed comment saying data was [n_samp-by-n_feat],
%                      added rand_init comment
%                      documented behavior when n_modes > numSamples
%                      changed rand_init selection when first selection
%                      fails to be unique.
%  Apr, 2008         Increased minimum initial stdev to 10*min_std 
%  May, 2009         implemeted fractional min_std
% May 2017  : updated min_std calculation
% Aug 2018  : changed warnings
%
%  Dr. Paul M Baggenstoss
%  NUWC, Newport, RI
%  p.m.baggenstoss@ieee.org
%

[n_feats,n_samples] = size(data);
if(nargin < 5), rand_init=1; end;
if nargin < 4 ,  min_std=.1; end;
min_std=min_std(:);
if nargin < 3 || isempty(names)
    for i = 1:n_feats, names{i} = ['Feat_' int2str(i)]; end
end;


if n_modes > n_samples
    msg = ['The number of modes can''t be more than the number of ' ...
        'samples... reducing n_modes to ' int2str(n_samples) '.'];
    warning(msg);
    n_modes = n_samples;
end

if(length(min_std) ~= n_feats),
    if length(min_std)==1,
        s=std(data,0,2);
        min_std = s*min_std;
    else,
        error('length(min_std) ~= n_feats');
    end;
end;

parm = [];
for (i_feat = 1:n_feats),
  parm.features(i_feat).name = names{i_feat};
  parm.features(i_feat).min_std = min_std(i_feat);
end;

% Determine sample mean and std of input data
% CAS: actually, this case is good for n_samples=1, too.
data_means = mean(data,2);
data_std = std(data,0,2);  % zero specifies default divisor (n-1)

% test for ill-conditioning
% if( any( abs(data_means) > 10000*data_std) || ...
%             max(data_std) > 5000*min(data_std)),
%        fprintf('***** INIT_GMIX: Red Alert!! Warning: Possible ill-conditioning: ');
%        fprintf(' [data_means(:) data_std(:)]=\n');
%        [data_means(:) data_std(:)]
%        
%        badIdx = find(abs(data_means) > 1000*data_std);
%        if isempty(badIdx)
%            fac = max(data_std)/min(data_std);
%            fprintf('The largest std is %d times the smallest!\n', fac);
%        else
%            fprintf('Mean > 1000*Std for Feature:  %s!\n', names{badIdx});
%        end
%           
%        fprintf('Your features are shit. You should scale them or remove means *****\n');
% end;

% CAS: this case is good for n_samples=1, too.
xmin = min(data,[],2);
xmax = max(data,[],2);

% implement fractional min_std
std_est = (xmax-xmin)/4;
min_std2=min_std;
ii=find(min_std<0);
min_std2(ii) = min_std(ii) .* std_est(ii);


% The initial value of covariances
starting_std = max((xmax-xmin)/n_modes/4,10*min_std2(:));

if( rand_init )
    % generate n_modes random numbers in [1,n_samples]
    idx = sort(1 + floor(rand(1,n_modes) * n_samples));
    if any(~diff(idx))
        [junk, p] = sort(rand(1, n_samples));   % sure-fire way
        idx = p(1:n_modes);
    end  
else	% non-random mode mean initialization
    idx = 1:n_modes;
end
means = data(:,idx);

%  Create the (Cholesky of) covariances and store the means
weight = 1/n_modes;
for (i_mode = 1:n_modes),
  parm.modes(i_mode).cholesky_covar = diag(starting_std);
  parm.modes(i_mode).mean = means(:,i_mode);
  parm.modes(i_mode).weight = weight;
end;

%fprintf('init_gmix successful. Features:');
%for i_feat = 1:n_feats,
%   fprintf(' %s', parm.features(i_feat).name);
%end;
%fprintf('\n');

