function [ log_post, alldata_log_lkld, log_joint, log_lkld ] = compute_post_raylsherl( raylsherl, data, sanity_check )

% function [ log_post, alldata_log_lkld, log_joint, log_lkld ] = compute_post_raylsherl( raylsherl, data, sanity_check )
%
% Estimate posterior probabilities, joint probabilities and overall
% likelihood on "data", given the model "raylsherl".  If the optional
% flag "sanity_check" is set to 1 (default value), we check that all
% posterior probabilities sum to 1.
%
% This is based on the "Rayleigh + Shifted Erlang" model defined in:
%
% "Unsupervised Spectral Subtraction",
% by G. Lathoud, M. Magimai.-Doss, B. Mesot and H. Bourlard,
% to appear in Proceedings of ASRU 2005.
%
% The pdf of this model can be briefly written as:
% p( data | raylsherl ) = raylsherl.p0 * raylpdf( data, raylsherl.sigma )
%   + ( 1 - raylsherl.p0 ) * erlangpdf( data - raylsherl.sigma * raylsherl.erlang_t, raylsherl.lambda )
%
%
% INPUT ARGUMENTS
%
% raylsherl: a structure with 'p0', 'sigma', 'erlang_t', 'lambda' fields.
% raylsherl.p0        is the prior probability of silence.
% raylsherl.sigma     is the Rayleigh parameter.
% raylsherl.erlang_t  defines the shift of the Shifted Erlang.
% raylsherl.lambda    is the parameter of the Shifted Erlang.
%
% data: a 2-D matrix of real values (e.g. magnitude values in a spectrogram).
%
% sanity_check (optional): a binary flag (0 or 1). Default value is 1.
%
%
% OUTPUT ARGUMENTS
%
% log_post:   a 3D matrix of size [ size( data ) 2 ].
% log_post( :,:,1 ): estimated posterior probabilities of silence.
% log_post( :,:,2 ): estimated posterior probabilities of activity.
%
% alldata_log_lkld:  a scalar, the log likelihood of the whole data.
% ( product of all lklds).
%
% log_joint:  a 3D matrix of size [ size( data ) 2 ].
% log_joint( :,:,1 ): the joint probabilities of data and silence.
% log_joint( :,:,2 ): the joint probabilities of data and activity.
%
% log_lkld:   a 2D matrix of size [ size( data ) ],
%   the likelihood of the data given the model.
%
% 
% SEE ALSO:  fit_raylsherl
%
%
% -- by G. Lathoud, 2005-08-26
  
  if nargin < 2
    error( [ mfilename ' needs 2 input arguments.' ] );
  end
  
  if ~exist( 'sanity_check', 'var' )
    sanity_check = 1;
  end
  
  % Check that the model has all necessary parameters
  
  check_param( { 'p0', 'sigma', 'erlang_t', 'lambda' }, fieldnames( raylsherl ) );
  
  % Check the data
  
  size_data = size( data );

  if max( find( size_data > 1 ) ) > 2
    error( [ mfilename ': data must be a 2-D matrix at most.' ] );
  end
  
  if ~all( data(:) >= 0 )
    error( [ mfilename ': all data must be >= 0' ] );
  end
  
  z_ind  = find( data == 0 );
  nz_ind = find( data > 0 );
  
  % Prepare memory for output
  
  log_post  = repmat( -Inf, [ size_data( 1:2 ), 2 ] );
  log_joint = repmat( -Inf, [ size_data( 1:2 ), 2 ] );
  
  % ( 1 ) Estimate joint probabilities
  
  % Special case for zero data -> arbitrary value to activate the
  % silence distribution This way we avoid unnecessarily complicate
  % modelling of zeroes with a Dirac.
  %
  % For an example of complete, rigorous modelling, see the
  % "compute_post_2d2g1sr.m"
    
  tmp = log_joint( :,:,1 );
  tmp( z_ind ) = log( raylsherl.p0 ) + 0;
  log_joint( :,:,1 ) = tmp;
  
  % General case: >0 data
  
  % Rayleigh: p( data, silence )
  
  tmp = log_joint( :,:,1 );
  tmp( nz_ind ) = log( raylsherl.p0 ) + log( data( nz_ind ) ) - 2 * log( raylsherl.sigma ) - ( data( nz_ind ) .^ 2 ) / ( 2 * raylsherl.sigma ^ 2 );
  log_joint( :,:,1 ) = tmp;
  
  % Shifted Erlang: p( data, activity )

  grt_ind = find( data > raylsherl.erlang_t );
  if ~isempty( grt_ind )
    tmp = log_joint( :,:,2 );
    data_grt = data( grt_ind ) - raylsherl.erlang_t;
    tmp( grt_ind ) = log( 1 - raylsherl.p0 ) + 2 * log( raylsherl.lambda ) + log( data_grt ) - raylsherl.lambda * data_grt;
    log_joint( :,:,2 ) = tmp;
  end
  
  % ( 2 ) Estimate posteriors
  
  log_lkld = my_logsum_fast( log_joint( :,:,1 ), log_joint( :,:,2 ) );
  
  log_post = log_joint - repmat( log_lkld, [ 1 1 2 ] );
  
  if sanity_check
    % Sanity check
    if any( abs( my_logsum_fast( log_post(:,:,1), log_post(:,:,2 ) ) ) > 1e-10 )
      error( [ mfilename ' posteriors do not seem to sum to one!' ] );
    end
  end
  
  % ( 3 ) Lkld of the whole data
  % This is useful to detect convergence of EM training
  
  if nargin > 1
    
    alldata_log_lkld = mean( log_lkld(:) );
    
    % Sanity check
    if isnan( alldata_log_lkld )
      error( [ mfilename ': unexpected NaN for alldata_log_lkld.' ] );
    end
    
  end
  
  if sanity_check
    % Sanity check
    if any( isnan( log_post(:) ) )
      error( [ mfilename ': unexpected NaN(s) in log_post.' ] );
    end
    if any( +Inf == log_post(:) )
      error( [ mfilename ': unexpected +Inf(s) in log_post.' ] );
    end
  end
  