function init_mvn_model = determinit_gmm( M, data )

% Deterministic initialization of a GMM
  
  if nargin < 2
    error( [ mfilename ' needs 2 input arguments.' ] );
  end
  
  mu = mean( data, 2 );
  Sigma = cov( data.' );
  
  % Make sure Sigma is well conditionned

  [ R, p ] = chol( Sigma );

  if p ~= 0
    Sigma = diag( diag( Sigma ) );
    [ R, p ] = chol( Sigma );
  end


  if ( ( numel( Sigma ) == 1 ) & ( numel( mu ) > 1 ) ) | ( size( data, 2 ) < 2 * M ) | ( p ~= 0 )
  
    % Problematic case (usually not enough data)
  
    Sigma = diag( ones( size( mu ) ) );
    
    disp( [ mfilename ' could not do proper GMM initialization. Fell back on 1-Gaussian component.' ] );

    init_mvn_model = [];
    init_mvn_model.w = 1 ;
    init_mvn_model.mu = mu;
    init_mvn_model.Sigma = Sigma;
    init_mvn_model.diag_Sigma = 1;

    return;

  end

  % Usual case (enough data)

  Rinv = inv( R );

  % Find the point that is furthest from the mean
  % ( most likely the component associated with this point will disappear below,
  % so we temporarily increase M )

  M2 = 2 * M;
  
  point_ind = [];
  
  d_mu = sum( ( ( data - repmat( mu, 1, size( data, 2 ) ) ).' * Rinv ) .^ 2, 2 ).';
  [ tmp, point_ind( 1 ) ] = max( d_mu );
  
  for m = 2:M2

    % Find the point that is furthest from the mean and all previously selected points

    d = d_mu;
    
    for p = 1:m-1
      
      d = min( d, sum( ( ( data - repmat( data( :, point_ind( p ) ), 1, size( data, 2 ) ) ).' * Rinv ) .^ 2, 2 ).' );
      
    end

    [ tmp, point_ind( m ) ] = max( d );

  end

  % Split the data between the M2 selected points (similar to centroids in K-means)

  d = repmat( NaN, M2, size( data, 2 ) );
  
  for m = 1:M2
    
    d( m,: ) = sum( ( ( data - repmat( data( :, point_ind( m ) ), 1, size( data, 2 ) ) ).' * Rinv ) .^ 2, 2 ).';

  end

  [ tmp, m_ind ] = min( d, [], 1 );
  
  % Keep only the M most represented components
  
  [ hw, hx ] = hist( m_ind, 1:M2 );
  
  [ tmp, ind ] = sort( -hw );
  
  point_ind = point_ind( ind( 1:M ) );
  
  %%% et rebelote

  % Split the data between the M selected points (similar to centroids in K-means)

  d = repmat( NaN, M, size( data, 2 ) );
  
  for m = 1:M
    
    d( m,: ) = sum( ( ( data - repmat( data( :, point_ind( m ) ), 1, size( data, 2 ) ) ).' * Rinv ) .^ 2, 2 ).';

  end

  [ tmp, m_ind ] = min( d, [], 1 );

  %%% Initialize the various components

  clear init_mvn_model

  m_new = 0;

  for m = 1:M

    ind = find( m_ind == m );

    if numel( ind ) < 2

      % Drop this component      

    else

      % Mean and covariance matrix of the selected points

      m_new = m_new + 1;

      init_mvn_model( m_new ).w  = numel( ind );  % We'll normalize below
      init_mvn_model( m_new ).mu = mean( data( :, ind ), 2 );
      init_mvn_model( m_new ).Sigma = cov( data( :, ind ).' );
      
      [ R, p ] = chol( init_mvn_model( m_new ).Sigma );
      
      if p == 0
	init_mvn_model( m_new ).diag_Sigma = 0;
      else
	init_mvn_model( m_new ).diag_Sigma = 1;
	init_mvn_model( m_new ).Sigma = diag( diag( init_mvn_model( m_new ).Sigma ) );
      end
      
    end

  end

  % We may have less components than required, at this point

  M = numel( init_mvn_model );

  % Weight normalization

  w = [ init_mvn_model.w ];
  w = w / sum( w );

  for m = 1:M
    init_mvn_model( m ).w = w( m );
  end
  