function y = log_cumsum( x )
  
  if nargin < 1
    error( [ mfilename ' needs 1 input argument.' ] );
  end
  
  size_x = size( x );
  
  % Flatten everything to 2-D
  % ( we'll cumsum along dimension one only )
  x = reshape( x, [ size_x( 1 )  prod( size_x( 2:end ) ) ] );
  
  y = x;
  
  for a = 2:size_x( 1 )
    
    ind = find( x( a,: ) ~= -Inf );
    
    if ~isempty( ind )
      y( a, ind ) = x( a, ind ) + log( 1 + exp( y( a-1, ind ) - x( a, ind ) ) );
    end
    
    ind = find( x( a, : ) == -Inf );
    y( a, ind ) = y( a-1, ind );
    
  end
  
  
  
  % Back to original dimensionality
  y = reshape( y, size_x );