function BF = BayesFactor(S,F_cell,D_cell,delta_vec,h_tilde,alf),
% BF = BayesFactor(S,F_cell,D_cell,delta_vec,h_tilde,alf),
%
% This code evaluates the Bayes factor U_{BF} in eq. (24) when the data
% is real and eq. (38) when it is complex. 
%
% Input
%   S: (p_rx x p_rx), sample convariance matrix.
%   p_rx: scalar, size of S.
%   F_cell: cell array, H0 parameters.
%   D_cell: cell array, H0 parameters.
%   delta_vec: (p_rx x 1), Gamma distribution parameters.
%   h_tilde: scalar, H1 parameters.
%   alf: scalar, H1 parameters, for the denominator in the Bayes Factor
%
% Output
%   BF: scalar, Bayes Factor as defined in eq. (24).
%
% Reference:
% M. H. Al-Ali and K.C. Ho, "Objective Bayesian approach for binary hypothesis
% testing of multivariate Gaussian observations," IEEE Trans. Inf. Theory,
% vol. 69, no. 2, pp. 1337-1354, Feb. 2023.
%
% M. H. Al-Ali and K. C. Ho, Dec. 12. 2022
%
%       Copyright (C) 2022
%       Compputer Engineering Department
%       University of Basrah
%       Basrah Province, Iraq
%       mohannd.khalaf@uobasrah.edu.iq
%
%       Copyright (C) 2022
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

ComplexData=(sum(sum(abs(imag(S))))>0);
if (ComplexData)
    Ftr=1;
else
    Ftr=0.5;
end;

% size of sample covariance matrix
p_rx = size(S,1);

% Calculate the Bayes factor Eq. (24)
S_cell = cell(p_rx, 1);
S_cell_tilde = cell(p_rx - 1, 1);
S_cell{1} = S(1,1);
c = zeros(p_rx, 1);
eta_tilda = zeros(p_rx, 1);
eta_tilda(1) = S(1,1);

for k = 2: p_rx
    s_vec = S(1:k-1,k);
    S_inv = S_cell{k - 1}^(-1);
    S_cell{k} = [S_cell{k-1} s_vec; s_vec' S(k,k)];
    S_cell_tilde{k-1} = S_cell{k-1} + (F_cell{k}*D_cell{k}*F_cell{k})^(-1);
    c(k) = s_vec'*(S_inv - S_cell_tilde{k-1}^(-1))*s_vec;
    eta_tilda(k) = S(k,k) - s_vec'*S_inv*s_vec;
end

eta_ratio = 1;
alf_max = max(max(alf), h_tilde);   % Used to scale exponents of the detector

for m = 1: p_rx
    num_var = (Ftr*(eta_tilda(m) + c(m)) + delta_vec(m))^(h_tilde/alf_max);
    den_var = (Ftr*eta_tilda(m))^(alf(m)/alf_max);
    eta_ratio = eta_ratio*(num_var/den_var);
end

BF = eta_ratio;
