%% Multivariate analysis
% Vine copula modelling

clear; 
clc;
close all;

addpath 'C:\Users\Menno\OneDrive\Documenten\TU Delft Civiele Techniek\Thesis\Coding\matlab scripts\waitbarParfor';

%% loading the data
datafit = readtable('Offshore GenoaExtremeWaves.txt');

Hs=datafit{:,2};
Tp=datafit{:,3};
Dir=datafit{:,4};
w_spd = datafit{:,5};
w_dir = datafit{:,6};

Hs_dens = ksdensity(Hs,Hs,'function','cdf');
Tp_dens = ksdensity(Tp,Tp,'function','cdf');
Dir_dens = ksdensity(Dir,Dir,'function','cdf');
w_spd_dens = ksdensity(w_spd,w_spd,'function','cdf');
w_dir_dens = ksdensity(w_dir,w_dir,'function','cdf');

datafit = [Hs_dens Tp_dens Dir_dens w_spd_dens w_dir_dens];

%% tree-equivalent vines
D = [    1 1 2 3 1; % V6
         0 2 1 2 2;
         0 0 3 1 3;
         0 0 0 4 4;
         0 0 0 0 5
         ]; 
B1 = [   1 1 1 1 2; % V7
         0 2 2 2 1;
         0 0 3 3 3;
         0 0 0 4 4;
         0 0 0 0 5
         ];
B2 = [   1 1 2 1 1; % V8
         0 2 1 2 2;
         0 0 3 3 3;
         0 0 0 4 4;
         0 0 0 0 5
         ];
B3 = [   1 1 1 1 2; % V8
         0 2 2 3 1;
         0 0 3 2 3;
         0 0 0 4 4;
         0 0 0 0 5
         ];
B0 = [   1 1 1 1 1; % V9
         0 2 2 2 3;
         0 0 3 3 2;
         0 0 0 4 4;
         0 0 0 0 5
         ];  
C = [    1 1 1 1 1; % V10
         0 2 2 2 2;
         0 0 3 3 3;
         0 0 0 4 4;
         0 0 0 0 5
         ];
  
     
%% permute every vine
EqC = cat(3,D,B1,B2,B3,B0,C);               % 3D array of all equivalence classes
n_EqC = size(EqC,3);                        % number of equivalence classes
datafit = [Hs_dens Tp_dens Dir_dens w_spd_dens w_dir_dens];          % data to be fitted
permutation = perms(1:1:size(datafit,2));   % all possible permutations of 5 nodes

% fit all possible vines
ppm = ParforProgressbar(length(permutation));
parfor ii = 1:length(permutation) % loop for every permutation
    % vary the permutation of the data each iteration
    data_perm(:,:,ii) = [datafit(:,permutation(ii,1)) datafit(:,permutation(ii,2)) datafit(:,permutation(ii,3)) datafit(:,permutation(ii,4)) datafit(:,permutation(ii,5))];
    
    for jj = 1:n_EqC % loop for every equivalence class
        [theta,~,sll] = ssp(data_perm(:,:,ii),EqC(:,:,jj),{'AIC'});
        AIC(jj,ii) = aicbic(sll,length(nonzeros(cell2mat(theta(1:16)))));
    end
    ppm.increment();
end
delete(ppm);

%% boxplot of tree equivalent vines
V6_AIC = uniquetol(AIC(1,:),0.0000001);
V7_AIC = uniquetol(AIC(2,:),0.0000001);
V8_AIC = [ uniquetol(AIC(3,:),0.0000001) uniquetol(AIC(4,:),0.0000001) ];
V9_AIC = uniquetol(AIC(5,:),0.0000001);
V10_AIC = uniquetol(AIC(6,:),0.0000001);

AIC_plot = [ V6_AIC V7_AIC V8_AIC V9_AIC V10_AIC];

savepath = 'C:\Users\Menno\OneDrive\Documenten\TU Delft Civiele Techniek\Thesis\Images\4-MultivariateModelling\';
grp = [ zeros(1,length(V6_AIC)), ones(1,length(V7_AIC)),2*ones(1,length(V8_AIC)),3*ones(1,length(V9_AIC)),4*ones(1,length(V10_AIC)) ];

figure
boxplot(AIC_plot,grp,'Labels',{'V6','V7','V8','V9','V10'})
ylabel('AIC')
figname = ['VineBoxplot.png'];
saveas(gcf,fullfile(savepath, figname)) 

%%
% AIC_list = zeros(5,180);
AIC_list(1,1:length(V6_AIC)) = V6_AIC;
AIC_list(2,1:length(V7_AIC)) = V7_AIC;
AIC_list(3,1:length(V8_AIC)) = V8_AIC;
AIC_list(4,1:length(V9_AIC)) = V9_AIC;
AIC_list(5,1:length(V10_AIC)) = V10_AIC;

q=uniquetol(AIC,0.0000001);

for ii = 1:480
    [row(ii),~] = find(AIC_list==q(ii));
    row(ii) = row(ii) + 5;
end

VinePerformance = [round(q,1) row'];

%% best fit
minAIC = min(AIC(:));
[row,col] = find(AIC==minAIC);


toc