function y = my_logsum_fast( log_x, opt_log_y )

% function y = my_logsum_fast( log_x )
%
% log_x : vector ( or matrix )containing log(xi)
% y     : scalar ( or row ), estimation of log(sum over i of xi)
%
% In case of a matrix we sum each column.
%
% Note that :
%  my_logsum_fast( [] ) returns -Inf
%  my_logsum_fast( [ -Inf -Inf ] ) returns -Inf
%  my_logsum_fast( [ -Inf -Inf 2 ] ) returns 2
%
% An explanation of the implementation is given in "logsum_matlab.pdf".
%
%
% function y = my_logsum_fast( log_x, log_y )
%
% Term-by-term log summation of two matrices of same dimension and
% size. Returns Y, a matrix of same size as LOG_X and LOG_Y:
% y = log( exp( log_x ) + exp( log_y ) )
%
%
% -- by G. Lathoud, 2005-08-26
  
  if nargin < 1
    error( 'logsum needs at least one input argument' );
  end

  if any( isnan( log_x(:) ) )
    error( [ mfilename ' log_x must not contain any NaN value!' ] );
  end
  
  if nargin > 1
    
    if any( isnan( opt_log_y(:) ) )
      error( [ mfilename ' opt_log_y must not contain any NaN value!' ] );
    end
    
    % Two matrices
    siz_x = size( log_x );
    siz_y = size( opt_log_y );
    if length( siz_x ) ~= length( siz_y )
      error( 'logsum: the two matrices must have same dimension and same size!' );
    end
    if ~all( siz_x == siz_y )
      error( 'logsum: the two matrices must have same dimension and same size!' );
    end

    if ~any( siz_x )
      % Empty matrices
      y = log_x;
      return;
    end
    
    M = max( log_x, opt_log_y );
    y = M + log( 1 + exp( min( log_x, opt_log_y ) - M ) );

    ind = find( ( log_x == -Inf ) & ( opt_log_y == -Inf ) );
    y( ind ) = -Inf;

    ind = find( ( log_x == +Inf ) & ( opt_log_y == +Inf ) );
    y( ind ) = +Inf;
    
    % Return output value
    return;
    
  end
  
  siz = size( log_x );
  
  n = length( find( siz > 1 ) );
  
  if n > 2
    error( [ mfilename ': can only deal with 2-D matrices.' ] );
  end
  
  if n == 0
    
    % Nothing to do

    if isempty( log_x )
  
      % Same convention as MATLAB sum( [] ) returning 0.
      y = -Inf;
    
    else
      
      y = log_x;
    
    end
    
    return;
  
  end
  
  if n == 1
    % Vector
    log_x = log_x(:);
  end
  
  % General case: matrix
  n_terms = size( log_x, 1 );
  
  % Start with the smallest value
  log_x = sort( log_x, 1 );
  
  % We do it in a hierarchical manner
  % At each iteration (k = 0...k_max),
  % we sum as many pairs of terms as possible
  k_max = log2( n_terms );

  % An explanation of the implementation is given in "logsum_matlab.pdf".
  
  for k = 0:k_max
    
    p_k = ceil( n_terms * (2^-k) );
    q_k = floor( p_k * 0.5 );
    t_k = 1 + ( 0:q_k-1 ) * 2 ^ (k+1);
    u_k = t_k + 2^k;

    tmp_t = log_x( t_k, : );
    tmp_u = log_x( u_k, : );
    
    ind1  = find( ( tmp_t == -Inf ) & ( tmp_u == -Inf ) );
    ind2  = find( ( tmp_t == +Inf ) & ( tmp_u == +Inf ) );
    ind3  = find( ( tmp_t == +Inf ) & ( tmp_u == -Inf ) );
    ind4  = find( ( tmp_t == -Inf ) & ( tmp_u == +Inf ) );
    
    tmp_t = tmp_u + log( 1 + exp( tmp_t - tmp_u ) );
    
    tmp_t( ind1 ) = -Inf;
    tmp_t( ind2 ) = +Inf;
    tmp_t( ind3 ) = +Inf;
    tmp_t( ind4 ) = +Inf;
    
    log_x( t_k, : ) = tmp_t;
    
  end
  
  % Return the result
  y = log_x( 1,: );