% Script to test Figure 2
% Comment/uncomment the relevant lines below depending on which
% regularization implementation you would like to test.


MC = input('Enter the number of MC trials to be used: ');
useSITB = 1; % 
%useSITB = 0; % Uncomment as desired
if useSITB
    disp('SITB will be used for the Regularized estimates')
else
    disp('The special code will be used for the Regularized estimates')
end
MaxOrder = 30;
was = warning;
warning('off')
for j = 1:MC
    % Generation of data:
    disp(['Generate data, Set number: ',int2str(j)])
    [Data,TestData,Syst] = DataGeneration_Sect7_3;
    LLData{j}=Data; LLTestData{j}=TestData; LLSyst{j}=Syst; % This is to store
        % the generated data and systems
    disp('Estimate Pem-oracle and Pem-AIC')
    [m_o,m_a,m_ac,m_b] = armax_with_oracle_and_with_AIC(Data,TestData,MaxOrder);
    disp('Estimate Regularized ARX')
    if useSITB
    m_DC = arx_with_regularization_with_SITB(Data,[40 40 1],'DC');
    m_TC = arx_with_regularization_with_SITB(Data,[40 40 1],'TC');
    m_SS = arx_with_regularization_with_SITB(Data,[40 40 1],'SS');
    else
    m_DC = arx_with_regularization_special_code(Data,[50 50 1],'DC');
    m_TC = arx_with_regularization_special_code(Data,[50 50 1],'TC');
    m_SS = arx_with_regularization_special_code(Data,[50 50 1],'SS'); 
    end
    
    disp('Estimate Pem-CV')
    m_CV = armax_with_CV(Data,MaxOrder);
    
    for ph=1:20
        [~,fitt] = compare(TestData,m_o,m_TC,m_SS,m_DC,m_CV,m_b,m_a,m_ac,ph,'ini','z');
        fit(j,ph,1:8)=max(cell2mat(fitt),-inf);
    end
end
disp('Show Results:')
fitph=squeeze(mean(fit,1));
figure(1) %Figure 2 left
plot([1:20],fitph(:,1),'r.',[1:20],fitph(:,2),'b*',[1:20],fitph(:,3),'bo',...
    [1:20],fitph(:,4),'bd',[1:20],fitph(:,5),'r+',[1:20],fitph(:,6),'m>',[1:20],fitph(:,8),'b^')
axis([0 20 -40 100]),grid
legend('PEM-Oracle','TC','SS','DC','PEM+CV','PEM+BIC','PEM+AICc'),shg 
figure(2) % Figure 2 right
text=['PEM-Or';'TC    ';'SS    ';'DC    ';'PEM-CV';'PEM+BI';'PEM+A ';'PEM-Ac'];

boxplot(squeeze(fit(:,1,[1:8])),text),grid, axis([0.5 6.5 0 100]),shg
warning(was);