% Script to test the Example in Sect 7.2
% Comment/uncomment the relevant lines below depending on which
% regularization implementation you would like to test.

useSITB = 1; % 
%useSITB = 0; % Uncomment as desired
MC = input('Enter the number of MC trials to be used: '); 
MaxOrder=30; % could be changed
if useSITB
    disp('SITB will be used for the Regularized estimates')
else
    disp('The special code will be used for the Regularized estimates')
end
 was=warning;
 warning('off')
for j = 1:MC
    % Generation of data:
    disp(['Generate data, Set number: ',int2str(j)])
    [EstData,TestData_w,TestData_u,Syst] = DataGeneration_Sect7_2;
    EstData=[iddata(zeros(200,1),zeros(200,1));EstData]; % This is to enforce 
       %the zero initial conditions for the regularized estimated
    LLSyst{j}=Syst; LLED{j}=EstData; LLTDw{j}=TestData_w; LLTDu{j}=TestData_u;
     % The above is to store the generated data and systems
   % Estimate models
    
    warning('off')
    for dlength=[1:2]% 1 is the short data record and 2 is the long one
        if dlength == 1
            Data = EstData(1:400);
        else
            Data = EstData;
        end
        disp('Estimate OE-oracle and OE-CV')
        
        mo_CV{dlength} = oe_with_CV(Data,MaxOrder);
        [mo_Or1{dlength},mo_Or2{dlength}] = ...
            oe_with_oracle_and_with_AIC(Data,TestData_w,TestData_u,MaxOrder);
         
        disp('Estimate Regularized FIR')
        if useSITB
            m_DC{dlength} = arx_with_regularization_with_SITB(Data,[0 170 1],'DC');
            m_SS{dlength} = arx_with_regularization_with_SITB(Data,[0 170 1],'SS');
            m_TC{dlength} = arx_with_regularization_with_SITB(Data,[0 170 1],'TC');
        else
            
            m_DC{dlength} = arx_with_regularization_special_code(Data,[0 200 1],'DC');
            m_SS{dlength} = arx_with_regularization_special_code(Data,[0 200 1],'SS');
            m_TC{dlength} = arx_with_regularization_special_code(Data,[0 200 1],'TC');
        end
    end
    % Evaluate the models using the test data
    % Think on the four corners of the subplots in Figure 1m when
    % interpreting fit{k,i}: 
    % {1,1} is the upper left, {1,2} is the upper right, {2,1} is the
    % lower left and {2,2} is the lower right
    for i=1:2 % Short Long
        for k=1:2 %1:white input test data, 2: system input testdata
            if k==1
                TestData = TestData_w;
            else
                TestData = TestData_u;
             end
            [~,fitt{k,i}] = compare(TestData,mo_Or1{i},m_TC{i},m_SS{i},...
                m_DC{i},mo_CV{i},mo_Or2{i},'ini','z');
            fit{k,i}(j,1:6)=max(cell2mat(fitt{k,i}),-inf);
        end 
    end
end
disp('Show Results:')
text=['OE+Or1';'TC    ';'SS    ';'DC    ';'OE+CV ';'OE+Or2'];
subplot(2,2,1)
boxplot(fit{1,1},text)
grid
title('short data; white input test')
axis([0.5 6.5 -20 100])
subplot(2,2,2)
boxplot(fit{1,2},text)
grid
title('Long data; white input test')
axis([0.5 6.5 -20 100])
subplot(2,2,3)
boxplot(fit{2,1},text),grid
title('Short data: BL input test')
axis([0.5 6.5 40 100])
subplot(2,2,4)
boxplot(fit{2,2},text),grid
title('Long Data: BL input test')
axis([0.5 6.5 40 100]),shg
MEAN_FITS=[mean(fit{1,1});mean(fit{1,2});mean(fit{2,1});mean(fit{2,2})];

% The Table below the plot:
disp('    OE+Or1     TC       SS         DC       OE+CV    OE+Or2')   

MEAN_FITS
disp('Table 1: Rows SD+WT; LD+WT; SD+BLT; LD+BLT')
warning(was);
 