function  [pd_ana,avg_prob_det]  = PdAnalyticalAndSimulation(m_tx, p_rx, K, lp_class, Prob_fa_vec, snr_db, rho, pow_gen)
% [pd_ana,avg_prob_det]  = PdAnalyticalAndSimulation(m_tx, p_rx, K, lp_class, Prob_fa_vec, snr_db, rho, pow_gen)
%
% This code computes ROC the analytic and simulation ROC curves for the 
% approximated U_{BF} detector in eq. (28).
%
% Input
%   m_tx: scalar, number of Tx antennas.
%   p_rx: scalar, number of Rx antennas.
%   K: scalar, number of received samples.
%   lp_class: scalar, index for the objective prior, possible values are
%       1,..,5.
%   Prob_fa_vec: (P x 1), P values of probability of false alarm.
%   snr_db: scalar, SNR value in dB.
%   rho: scalar, correlation parameter to produce eq. (40) in the reference paper. 
%   pow_gen: (p_rx x 1), Diagonal values of Omega_0 in the log scale.
%
% Output
%   pd_ana: (P x 1), Analytical probability of detection at the prob. of false alarms.
%   avg_prob_det: (P x 1), Simulated probability of detection at the prob. false alarms.
%
% 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
%
 

% -- Objective prior parameters --
eps_val = .01;                      % A small parameter to satisfy eq. (27)
m = 1:p_rx;

if lp_class == 1                                
    a_m = p_rx - m;                 % Jeffreys Prior
    h = 1/2 + eps_val;
elseif lp_class == 2
    a_m = p_rx - m + 1;             % Independence jeffreys Prior
    h = 0 + eps_val;
elseif lp_class == 3
    a_m = 2 - m;                    % Geisser and Cornfield Prior
    h = (p_rx - 1)/2 + eps_val;
elseif lp_class == 4
    a_m = m;                        % Right Haar Measure
    h = 0 + eps_val;
elseif lp_class == 5
    a_m = ones(1, p_rx);            % Reference Prior
    h = 0 + eps_val;
end

% -- Simulation parameters --
num_rlzn = 5e5;                      % Number of realizations
fa_length = length(Prob_fa_vec);
snr_linear = 10^(snr_db/10);

% -- Covariance matrix for null hypothesis --
Cov_n = zeros(p_rx, p_rx);

for i = 1: p_rx
    for j = 1: p_rx
        if i <= j
            Cov_n(i, j) = rho^(j - i);
        else
            Cov_n(i, j) = rho^( i - j);
        end
    end
end

% -- noise power --
pow_noisedB = pow_gen(1: p_rx);
pow_noise = 10.^(pow_noisedB/10);
D = diag(pow_noise);
Omega_0 = (D.^(.5)*Cov_n*D.^(.5));

H = randn(p_rx, m_tx);          % Generate one Ch matrix for all the realizations

Cov_H = H*H';
normalize_cov_H = Cov_H/trace(Cov_H);
sig_pow = snr_linear*trace(Omega_0);
Omega_1 = sig_pow*normalize_cov_H + Omega_0;
Omega0_cell = cell(p_rx,1);
Omega1_cell = cell(p_rx,1);
Omega0_cell{1} =  Omega_0(1, 1);
Omega1_cell{1} =  Omega_1(1, 1);

for k = 2: p_rx
    w0_vec = Omega_0(1:k-1,k);
    w1_vec = Omega_1(1:k-1,k);
    Omega0_cell{k} = [Omega0_cell{k-1} w0_vec; w0_vec' Omega_0(k,k)];
    Omega1_cell{k} = [Omega1_cell{k-1} w1_vec; w1_vec' Omega_1(k,k)];
end

% ----- Analytical Pd evaluation -----

% -- Moment under null model eq. (34) -- 
Th_val_ana_B = zeros(fa_length, 1);
moment_fa = zeros(3, 1);
moment_pd = zeros(3, 1);

for q = 1:3
    
    coeff_fa = 1;
    coeff_pd = 1;
    
    for i = 1: p_rx
        a = (h + (a_m(i) - 1)/2);
        
        if i == 1
            coeff_fa = coeff_fa*(2*det(Omega0_cell{i}))^(q*a)...
                *gamma(q*a)/beta(q*a, (K - i)/2);
            coeff_pd = coeff_pd*(2*det(Omega1_cell{i}))^(q*a)...
                *gamma(q*a)/beta(q*a, (K - i)/2);
        else
            coeff_fa = coeff_fa*(2*det(Omega0_cell{i})/det(Omega0_cell{i - 1}))^(q*a)...
                *gamma(q*a)/beta(q*a, (K - i)/2);
            coeff_pd = coeff_pd*(2*det(Omega1_cell{i})/det(Omega1_cell{i - 1}))^(q*a)...
                *gamma(q*a)/beta(q*a, (K - i)/2);
        end
        
    end
    
    moment_fa(q) = coeff_fa;
    moment_pd(q) = coeff_pd;
    
end

sigma_fa = sqrt(moment_fa(2) - moment_fa(1)^2);
Skew_fa = (moment_fa(3) - 3*moment_fa(1)*sigma_fa^2 - moment_fa(1)^3)/sigma_fa^3;
alpha_fa = 4/Skew_fa^2;
beta_fa = sqrt(alpha_fa)/sigma_fa;
s_fa = alpha_fa/beta_fa - moment_fa(1);

% -- Moments under signal presence model --
sigma_pd = sqrt(moment_pd(2) - moment_pd(1)^2);
Skew_pd = (moment_pd(3) - 3*moment_pd(1)*sigma_pd^2 - moment_pd(1)^3)/sigma_pd^3;
alpha_pd = 4/Skew_pd^2;
beta_pd = sqrt(alpha_pd)/sigma_pd;
s_pd = alpha_pd/beta_pd - moment_pd(1);
pd_ana = zeros(fa_length, 1);

for lp_fa = 1:fa_length
    
    Th_val_ana_B(lp_fa) = (1/beta_fa)*gammaincinv(Prob_fa_vec(lp_fa), alpha_fa,'upper') - s_fa;
    pd_ana(lp_fa) = gammainc(beta_pd*(Th_val_ana_B(lp_fa) + s_pd), alpha_pd, 'upper') ;
    
end

% ----- Simulation Pd Evaluation -----

% -- Simulated threshold --
Th_val_sim = zeros(fa_length, 1);

for lp_fa = 1: fa_length
    
    [Th_val_sim(lp_fa), ~, ~] = DtnThFullMtxApprox(p_rx, Omega_0, Prob_fa_vec(lp_fa), K, h, num_rlzn, lp_class);
    
end

Prob_det = zeros(num_rlzn,fa_length);

for lp_real = 1:num_rlzn

    % Generate model data
    Y = mvnrnd(zeros(1, p_rx), Omega_1, K);     % Generates an K-by-p_rx matrix of random vectors

    % Sample covariance matrix
    S = Y'*Y;

    S_cell = cell(p_rx,1);
    S_cell{1} = S(1,1);
    eta_tilda = zeros(p_rx,1);
    eta_tilda(1) = S(1,1);
   
    % Bayes factor calculation for right Haar measure
    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)];
        eta_tilda(k) = S(k,k) - s_vec'*S_inv*s_vec;
    end
    
    eta_ratio = 1;
    
    for m = 1: p_rx
        eta_ratio = eta_ratio*eta_tilda(m)^(h + (a_m(m) - 1)/2);
    end
    
    Bayes_factor = (eta_ratio);
    
    % Detection decision 
    for lp_fa = 1:fa_length
        Prob_det(lp_real, lp_fa) = Bayes_factor >= Th_val_sim(lp_fa);
    end

end

avg_prob_det = mean(Prob_det)';


% =================================================================================================
function [Th_val, Out_Pfa, BF] = DtnThFullMtxApprox(p_rx, Omega_0, Pfa, K, h, num_ensembles, lp_class) 
% [Th_val, Out_Pfa, BF] = DtnThFullMtxApprox(p_rx, Omega_0, Pfa, K, h, num_ensembles, lp_class) 
%
% The code provides the threshold value for the approximated U_{BF} test in
% eq. (28) at a given prob of false alarm, for rifht Haar measure objective 
% prior for H1.
%
% Input
%   p_rx: scalar, number of Rx antennas.
%   Omega_0: (p_rx x p_rx): Covariance matrix in the null hypothesis,
%   Pfa: scalar, desired prob of false alarm.
%   K: scalar, number of received samples.
%   h: parameter for the prior of H1.
%   num_ensembles: scalar, number of realizations.
%   lp_class: scalar, index for the objective prior, possible values are
%       1,..,5.
%
% Output
%   Th_val: scalar, threshold value for the given Pfa.
%   Out_Pfa: scalar, Simulated prob. of false alarm to match with the desired Pfa.
%   BF ( num_ensembles x 1):  Bayes factor values.

% 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
%

% -- Objective prior parameters --
m = 1:p_rx;
if lp_class == 1                                                          % Possible values for a_m: su_rx-k, su_rx+1-k, 2-k, k, or 1. 
       a_m = p_rx - m;                                                    % Jeffreys Prior.
elseif lp_class == 2
       a_m = p_rx - m + 1;                                                % Independence jeffreys Prior.
elseif lp_class == 3
       a_m = 2 - m;                                                       % Geisser and Cornfield Prior.
elseif lp_class == 4 
       a_m = m;                                                           % Right Haar Measure.
elseif lp_class == 5
       a_m = ones(1, p_rx);                                               % Ref Prior.
end

% -- Null model generation --
BF = zeros(num_ensembles, 1);

for lp_real = 1:num_ensembles
          
           Y = mvnrnd(zeros(1, p_rx), Omega_0, K);                 % Generates an K-by-su_rx matrix of random vectors 
                                                                    % chosen from the multivariate normal distribution, 
           % Sample covariance matrix
           S = Y'*Y;

           % Bayes factor calculation
           S_cell = cell(p_rx, 1);
           S_cell{1} = S(1,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)];
               eta_tilda(k) = S(k,k) - s_vec'*S_inv*s_vec;
           end
           
               eta_ratio = 1;
               
               for m = 1: p_rx
                   eta_ratio = eta_ratio*(eta_tilda(m))^(h + (a_m(m) - 1)/2); 
               end   
               
               BF(lp_real, 1) = (eta_ratio);  
end

% -- Determining threshold --
nbins = 1e4;            % step resolution to determine threshold 
thresh_low = min(BF);
thresh_hi  = max(BF);
thresh_steps = linspace(thresh_low,thresh_hi,nbins);
pfa_count = zeros(1,nbins);

for k = 1:nbins
    thresh = thresh_steps(k);
    pfa_count(k) = sum(BF >= thresh);
end

pfa_sim = pfa_count/num_ensembles;
[err, indx] = min(abs(pfa_sim - Pfa));
Th_val = thresh_steps(indx);
Out_Pfa = pfa_sim(indx);
