function [ frame_log_post, sector_log_post, alldata_log_plkld, frame_log_plkld, frame_log_pjoint, sector_log_pjoint_cond ] = compute_post_2d2g1se( model, data, sanity_check )

% function [ frame_log_post, sector_log_post, alldata_log_plkld, frame_log_plkld, frame_log_pjoint, sector_log_pjoint_cond ] = compute_post_2d2g1se( model, data, [ sanity_check ] )
%
% ( sanity_check: default is 0.
%   If set to 1, then the script verifies that posteriors sum to one )
%
% Compute posteriors, partial likelihoods and partial joints for some
% data given a 2d2g1se model, which is defined as:
%
% p( x | sector_class, frame_class, model ) =
%   if frame_class = silence:
%     if sector_class = silence:
%       w00d * Dirac_0( x ) + w00g * GammaPdf( x, gamma00, beta00 )
%     if sector_class = activity:
%       0
%   if frame_class = activity: 
%     if sector_class = silence:
%       w01d * Dirac_0( x ) + w01g * GammaPdf( x, gamma01, beta01 )
%     if sector_class = activity:
%       if x <= gamma01 * beta01:  0
%       if x >  gamma01 * beta01:  RicePdf( x - gamma01 * beta01, sigma11, V11 )
%
% The priors of each class are:
%   a0 = P( frame_class = silence )
%   a1 = P( frame_class = activity )
%   b0 = P( sector_class = silence  | frame_class = activity )
%   b1 = P( sector_class = activity | frame_class = activity )
%
% Notes:  a0 + a1 = 1, b0 + b1 = 1  and  w0d + w0g = 1
% gamma00 > 1, beta00 > 0, gamma01 > 1 and beta01 > 0.
% The "gamma > 1" constraint ensures that GammaPdf( 0, gamma, beta ) = 0.
% 
% See the NIST handbook for a definition of the Gamma distribution
% http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
%
% See Mathworld for a definition of the Rice distribution
% http://mathworld.wolfram.com/RiceDistribution.html
%
% All outputs are defined in the log domain:
% Compute (exact) posteriors  -> frame_log_post, sector_log_post
% Compute partial likelihoods -> frame_log_plkld
% Compute partial lkld of whole data -> alldata_log_plkld = mean( frame_log_plkld(:) ) 
% Compute partial joints             -> sector_log_pjoint_cond, frame_log_pjoint
% "partial" means that all Diracs are replaced with Kroneckers.
%
% -> result: frame_log_post and sector_log_post are exact,
%            frame_log_plkld is partial: log( product of finite terms only )
%
% Note: full lklds/joints are not computed because = +Inf as soon as
% at least one data sample in "data" is equal to zero.
%
% alldata_log_plkld is useful e.g. to compare
% the fit of the 2d2g1se model for two different sets of parameters,
% because of the following exact equality:
%    
%    alldata_plkld( model_1 )   alldata_lkld( model_1 )
%    ------------------------ = -----------------------
%    alldata_plkld( model_2 )   alldata_lkld( model_2 )
%
% Dirac terms from the right hand side have been normalized out in the
% left hand side. Their associated weights remain in the product (w00d, w01d).
%
%
%
% INPUT:
%
% model is a structure with fields  
% { 'a0', 'b0', 'w00d', 'gamma00', 'beta00', 'w01d', 'gamma01', 'beta01', 'sigma11', 'V11' }.
% 
% ( Implicitly:  a1 = 1 - a0, b1 = 1 - b0, w00g = 1 - w00d and w01g = 1 - w01d. )
%
% data is a 2-D matrix of >=0 values.
%   rows: sectors of space.
%   columns: time frames.
%
% OUTPUT
% 
% frame_log_post  is a 3-D matrix of size = [ 1   size( data,2 )  2 ]
%
% sector_log_post is a 3-D matrix of size = [   size( data )      3 ]
% sector_log_post( :,:,1 ) contains the posterior of ( silent sector, silent frame ) for each data sample,
% sector_log_post( :,:,2 ) contains the posterior of ( silent sector, active frame ) for each data sample,
% sector_log_post( :,:,3 ) contains the posterior of ( active sector, active frame ) for each data sample.
% 
% alldata_log_plkld is a scalar
%
% frame_log_plkld  is a 3-D matrix of size = [ 1   size( data,2 )  1 ]
%
% frame_log_pjoint is a 3-D matrix of size = [ 1   size( data,2 )  2 ] 
%
% sector_log_pjoint_cond is a 3-D matrix of size = [ size( data ) 3 ]
% sector_log_pjoint_cond( u,v,1 ) contains:
%   p( data(u,v), silent sector (u,v) | model, silent frame v ) / product of all Dirac terms
% sector_log_pjoint_cond( u,v,2 ) contains:
%   p( data(u,v), silent sector (u,v) | model, active frame v ) / product of all Dirac terms
% sector_log_pjoint_cond( u,v,3 ) contains:
%   p( data(u,v), active sector (u,v) | model, active frame v ) / product of all Dirac terms
%
% Note: 
%   alldata_log_plkld = mean( frame_log_plkld(:) )
%   frame_log_plkld  = logsum along 3rd dimension of frame_log_pjoint
%   frame_log_pjoint( :,:,1 ) = log( a0 ) + sum( sector_log_pjoint_cond( :,:,1 ),1 )
%   frame_log_pjoint( :,:,2 ) = log( a1 ) + sum( logsum( log( b0 ) + sector_log_pjoint_cond( :,:,2 ), log( b1 ) + sector_log_pjoint_cond( :,:,3 ) ), 1 )
%   frame_log_post  = frame_log_pjoint - repmat( frame_log_plkld, [1 1 3] )
%
%
% SEE ALSO fit_2d2g1se, generate_2d2g1se
  
  
  if nargin < 2
    error( [ mfilename ' needs at least 2 input arguments.' ] );
  end
  
  if nargin < 3
    sanity_check = 0;
  end
  
  
  % Check the model
  % model: Make sure all parameters are defined
  
  model = check_2d2g1se_model( model );
  
    
  % data: make sure it is a 2-D matrix
  
  size_data = size( data );
  
  if max( find( size_data > 1 ) ) > 2
    error( [ mfilename ': "data" must be a 2-D matrix.' ] );
  end
  
  % Check all data is >= 0
  
  if any( data(:) < 0 )
    error( [ mfilename ': all data must be >=0.' ] );
  end
  
  z_ind  = find( data == 0 );
  nz_ind = find( data > 0 );
  
  % Prepare memory
  
  frame_log_post    = repmat( NaN, [ 1 size_data(2) 2 ] );
  frame_log_plkld   = repmat( NaN, [ 1 size_data(2) 1 ] );
  frame_log_pjoint  = repmat( NaN, [ 1 size_data(2) 2 ] );
  sector_log_post   = repmat( NaN, [ size_data 3 ] );
  sector_log_pjoint_cond = repmat( NaN, [ size_data 3 ] );
  
  % Compute partial joint of data and silent sector given silent frame
  % "1d1g" = mixture of a Dirac at zero and a Gamma pdf with gamma00 > 1
  
  a_mat = sector_log_pjoint_cond( :,:,1 );
  gamma00 = model.gamma00;
  beta00  = model.beta00;
  
  a_mat( z_ind )  = log( model.w00d );  % Dirac replaced with Kronecker
  a_mat( nz_ind ) = log( model.w00g ) + ( gamma00 - 1 ) .* log( data( nz_ind ) ) - data( nz_ind ) / beta00 - gammaln( gamma00 ) - gamma00 * log( beta00 );
  
  sector_log_pjoint_cond( :,:,1 ) = a_mat;
  
  % Compute partial joint of data and silent sector given active frame
  % "1d1g" = mixture of a Dirac at zero and a Gamma pdf with gamma01 > 1
  
  a_mat = sector_log_pjoint_cond( :,:,2 );
  gamma01 = model.gamma01;
  beta01  = model.beta01;
  
  a_mat( z_ind )  = log( model.w01d );  % Dirac replaced with Kronecker
  a_mat( nz_ind ) = log( model.w01g ) + ( gamma01 - 1 ) .* log( data( nz_ind ) ) - data( nz_ind ) / beta01 - gammaln( gamma01 ) - gamma01 * log( beta01 );
  
  sector_log_pjoint_cond( :,:,2 ) = a_mat + log( model.b0 );
  
  % Compute partial joint of data and active sector given active frame
  % "1se" = one shifted Erlang, shift = gamma01 * beta01
  
  a_mat = sector_log_pjoint_cond( :,:,3 );
  
  t11 = model.gamma01 * model.beta01;
  beta11 = model.beta11;
  
  a_mat( find( data <= t11 ) ) = -Inf;
  
  grt_ind = find( data > t11 );
  tmp = data( grt_ind ) - t11;
  a_mat( grt_ind ) = log( tmp ) - tmp / beta11 - 2 * log( beta11 );
  
  sector_log_pjoint_cond( :,:,3 ) = a_mat + log( model.b1 );
  
  % Sanity check
  
  if any( sector_log_pjoint_cond(:) == +Inf )
    error( [ mfilename ' is +Inf-insane.' ] );
  end
  
  if any( isnan( sector_log_pjoint_cond(:) ) )
    error( [ mfilename ' is NaN-insane.' ] );
  end
  
  % Conditional lkld: given active frame
  
  sector_log_plkld_g1 = my_logsum_fast( sector_log_pjoint_cond( :,:,2 ), sector_log_pjoint_cond( :,:,3 ) );
  
  % Frame-level joints
  
  frame_log_pjoint( :,:,1 ) = log( model.a0 ) + sum( sector_log_pjoint_cond( :,:,1 ), 1 );
  frame_log_pjoint( :,:,2 ) = log( model.a1 ) + sum( sector_log_plkld_g1, 1 );
  
  % Frame-level lkld
  
  frame_log_plkld   = my_logsum_fast( frame_log_pjoint( :,:,1 ), frame_log_pjoint( :,:,2 ) );
  
  % By-product: partial lkld of the entire data
  
  alldata_log_plkld = mean( frame_log_plkld(:) );
  
  % Frame-level posteriors
  
  frame_log_post  = frame_log_pjoint - repmat( frame_log_plkld, [ 1 1 2 ] );
  
  
  
  % Sector-level posteriors
  
  sector_log_post( :,:,1 ) = repmat( frame_log_post( :,:,1 ), [ size( sector_log_post, 1 ), 1, 1 ] );
  
  sector_log_post_g1 = sector_log_pjoint_cond( :,:,2:3 ) - repmat( sector_log_plkld_g1, [ 1 1 2 ] );
  
  sector_log_post( :,:,2:3 ) = sector_log_post_g1 + repmat( frame_log_post( :,:,2 ), [ size( sector_log_post, 1 ), 1, 2 ] );
  
  if sanity_check
    tmp = my_logsum_fast( reshape( sector_log_post, [], 3 ).' );
    
    if any( abs( tmp ) > 1e-10 )
      error( [ mfilename ': red alert! Posteriors do not sum to one!' ] );
    end
    
    disp( [ mfilename ' seems to be sane !' ] );
    
  end
  