% U               = gallery('qmult',k);
% Lambda          = diag(rand(k,1)*2 + 1);
% V               = gallery('qmult',k);
% B               = U*Lambda*V';
B = [0.3175   -0.9259   -0.4594    0.7967   -0.3507    0.6462    0.3316   -0.2903  1.2321;
     0.4813    0.1237    0.5299    1.2479    0.5493   -0.5761   -0.0520    0.2569  -0.3199;
     -1.1326   -0.4329    0.6435    0.3616   -0.4940   -0.3085    0.3345   -0.0493  -0.2063;
     0.4087   -0.7299    0.8755   -0.6444    0.7951    0.4749    0.8641    0.0717  0.0270;
     0.1631    0.7741   -0.0572    0.1253   -1.1738    0.0493    0.0499    1.0300 0.8116;
     0.3521   -1.0594   -0.2978    0.0582   -0.5542    0.2866   -0.5187    0.5630 -0.3212;
     -0.3513    1.0874    0.5375   -0.0708   -0.4345    1.1174   -0.5680   -0.1517  -0.0000;
     -0.9420    0.0717   -0.5846    0.2159    0.7213    1.0357    0.5871    0.6148    -0.1014;
     -1.0494   -0.5347    0.2006    0.2077    0.5169   -0.3007   -1.1130    0.5811  0.4759;];

noise_power         = 0.2;
n                   = 100000;
k                   = 9;
eps_vals            = 0.5:0.05:1;
num_trials          = 200;
num_repetitions     = 10;
itr                 = 1;
amari_scores        = zeros(length(eps_vals), num_repetitions);
independence_scores = zeros(length(eps_vals), num_repetitions);
uncorrected_independence_score = zeros(length(eps_vals), num_repetitions);
random_projections  = randn(2*num_trials,k);


for eps = eps_vals
    fprintf("Iteration : %d\n", itr);
    for r=1:num_repetitions
        [Sigma,B,X,z,e]            = generate_data_belkin(n,k,noise_power,B);
        B_prime                    = eps*B + (1-eps)*eye(k);
        amari_scores(itr,r)        = compute_amari_score(B_prime,B);
        independence_scores(itr,r) = Delta_Score_Total(X,B*B',B_prime,random_projections,num_trials);
        uncorrected_independence_score(itr,r) = Delta_Score_Uncorrected(X,B*B',B_prime,random_projections,num_trials);
    end
    fprintf("Amari distance = %.5f, Corrected Independence Score = %.5f, Uncorrected Independence Score = %.5f\n", mean(amari_scores(itr,:)), mean(independence_scores(itr,:)), mean(uncorrected_independence_score(itr,:)));
    itr = itr + 1;
end

mean_amari_score = mean(amari_scores, 2);
mean_ind_score   = mean(independence_scores, 2);
std_ind_score    = std(independence_scores, 0, 2);
mean_un_ind_score   = mean(uncorrected_independence_score, 2);
std_un_ind_score    = std(uncorrected_independence_score, 0, 2);

figure;
errorbar(mean_amari_score', mean_ind_score', std_ind_score,'LineWidth',3,'Marker','o','MarkerFaceColor','r'); hold on;
% errorbar(mean_amari_score', mean_un_ind_score', std_un_ind_score,'LineWidth',2,'Marker','o','MarkerFaceColor','r');
xlabel("Amari distance from true mixing matrix",'FontWeight','bold','FontSize', 30);
ylabel("Independence Score",'FontWeight','bold','FontSize', 30);
title("Correlation of Independence score with Amari distance",'FontWeight','bold','FontSize',30);
% legend("Corrected Independence Score", "Uncorrected Independence Score");
ax=gca;
ax.FontSize = 30;
ax.FontWeight = "bold";
grid on;