clear all; %#ok<CLALL>
close all;
clc;
rng(0);
cd(fileparts(matlab.desktop.editor.getActiveFilename));

%Setup inputs
BurnInPeriod = 4000;
QAPReps = 10000 * 100;
NumFactorsK = 2;

%Load network from Gibbs sampler
load(strcat('.\Output\RActual_SmallSim_H1_Draws6000.mat'));

%Get data after burn in period
GibbsRAtoRA   = EstimatedNetwork.RAtoRA_Firm(:,:,BurnInPeriod+1:end);
GibbsRIFtoRAF = EstimatedNetwork.RIFtoRAF(:,:,BurnInPeriod+1:end);

%Calculate medians across draws after burn in period
GibbsRAtoRAMed   = median(GibbsRAtoRA,3);
GibbsRIFtoRAFMed = median(GibbsRIFtoRAF,3);

clearvars EstimatedNetwork BurnInPeriod;

%Load OCMT network
OCMT_Full = readtable(strcat('.\Output\FullCon_RActual_SmallSimGFEVc_F2Op0.1O1FALSEOd11Od22OL1FL1H1.csv'),ReadRowNames=true);
OCMT_Full.Properties.VariableNames = OCMT_Full.Properties.RowNames;

NumSeriesN = (size(OCMT_Full,2)-NumFactorsK)/2;

OCMTRAtoRA   = table2array(    OCMT_Full(1:NumSeriesN,                            1:NumSeriesN));
OCMTRIFtoRAF = table2array([   OCMT_Full(1:NumSeriesN,                            NumSeriesN+1:2*NumSeriesN+NumFactorsK)   ;
                               OCMT_Full(2*NumSeriesN+1:2*NumSeriesN+NumFactorsK, NumSeriesN+1:2*NumSeriesN+NumFactorsK)   ]);


%Calculate correlations with median Gibbs network & QAP test values for each case
format long;
tic;

fprintf('\n');
disp("RAtoRA");
[QAPTestCorrAvg, QAPTestpValue] = qapBayesTest(OCMTRAtoRA, GibbsRAtoRA, QAPReps);
disp(corr(GibbsRAtoRAMed(:), OCMTRAtoRA(:)));
disp(QAPTestCorrAvg);
disp(QAPTestpValue);
toc;

fprintf('\n');
disp("Full RIFtoRAF");
[QAPTestCorrAvg, QAPTestpValue] = qapBayesTest(OCMTRIFtoRAF, GibbsRIFtoRAF, QAPReps);
disp(corr(GibbsRIFtoRAFMed(:), OCMTRIFtoRAF(:)));
disp(QAPTestCorrAvg);
disp(QAPTestpValue);
toc;

fprintf('\n');
disp("RIFtoRAF without factors to focus on inter-panel member connections");
[QAPTestCorrAvg, QAPTestpValue] = qapBayesTest(OCMTRIFtoRAF(1:NumSeriesN,1:NumSeriesN), GibbsRIFtoRAF(1:NumSeriesN,1:NumSeriesN,:), QAPReps);
disp(corr(reshape(GibbsRIFtoRAFMed(1:NumSeriesN,1:NumSeriesN),[NumSeriesN^2 1]), reshape(OCMTRIFtoRAF(1:NumSeriesN,1:NumSeriesN),[NumSeriesN^2 1])));
disp(QAPTestCorrAvg);
disp(QAPTestpValue);
toc;

disp("Done!");