clear all;

%%% FINAL VERSION: 17 MARCH 2021 %%%
%%% SENSITIVITY ANALYSES OF THEOREM 6.7 FOR THE CONDITIONS
%%% (A) EQUAL ATTENTIVENESS AND (B) CONSTANT COLUMNS OF ONE OF THE INTERGROUP
%%% MATRICES (A_{12} OR A_{21})

%A SENSITIVITY ANALYSES IS PERFORMED ON A POPULATION OF 4 AGENTS
%CONSISTING OF TWO GROUPS OF TWO AGENTS

%Define the number of agents within each group
m1 = 2;
m2 = 2;

%Interaction within group 1 is chosen such that it shows polarization
%asymptotically (1.2 + 0.99 > 1 + 1, apply Theorem 6.1)
P11 = [1.2 -0.2; 0.01 0.99];

%Create a diagonal matrix that normalizes the sum of each row of P11 to 1
P11_norm = inv(diag(P11 * ones(m1,1)));

%Adjust P11 such that P11 is a 1-equal attentiveness matrix
P11 = P11_norm * P11;

%Construct P12 such that it is a 1-equal attentiveness matrix and has
%constant columns. Sample the entries of the columns randomly between [-0.5,0.5]
P12_vec = (rand(1,m2))-0.5;

%Normalize the vector such that the entries sum up to 1
P12_norm = P12_vec / (P12_vec * ones(m2,1));

%Construct P12 where each column equals the entries of P12_norm
P12 = ones(m1,1) * P12_norm;

%Construct P22 similar as P11 with slightly different parameters
P22 = [0.9 0.1; -0.2 1.2];
P22_norm = inv(diag(P22 * ones(m2,1)));
P22 = P22_norm * P22;

%Construct P21 as an 1-equal attentiveness matrix and entries sampled
%uniformly from between [-0.5,0.5]
P21 = (rand(m2,m1)-0.5)*2;
P21_norm = inv(diag(P21 * ones(m1,1)));
P21 = P21_norm * P21;

%Choose a C that leads to polarization
C = [1 -0.2; -0.01 1.2];

%Construct A with C, P11, P12, P21 and P22
A = [C(1,1)*P11 C(1,2)*P12; C(2,1)*P21 C(2,2)*P22];

%Save the size of A
N = size(A,1);

%Iteration steps
steps = 60;

%Initiate a matrix of 0s to store the eigenvalues in
v_eigenvalues = zeros(N,steps);

%CHOOSE A SCENARIO TO RUN
%scenario 1: no noise
%scenario 2: P12 dynamic noise increase by 0.02 per iteration
%scenario 3: P12, P21, P11, P22 attentiveness noise by 0.005 per iteration

scenario = 3;

for i = 1:1:steps
 
    if scenario == 1
        disp('Scenario 1: no noise')   
        %Choose P12 as a 1-equal attentiveness matrix with constant columns
        %with entries uniformly sampled from [-0.5,0.5]
        P12_vec = rand(1,m2)-0.5;
        P12_norm = P12_vec / (P12_vec * ones(m2,1));
        P12 = ones(m1,1) * P12_norm;

        %Construct P21 as an 1-equal attentiveness matrix and entries
        %sampled uniformly from [-0.5,0.5]
        P21 = (rand(m2,m1))*2;
        P21_norm = inv(diag(P21 * ones(m1,1)));
        P21 = P21_norm * P21;

        %Construct A
        A = [C(1,1)*P11 C(1,2)*P12; C(2,1)*P21 C(2,2)*P22]; 
    
        %SaveName = ['A_SensANoNoiseIteration' num2str(i)];
        %save(SaveName,'A');
        
        ExportDPI = 300;
        ExportName = 'SensANoNoise300.png';
        
    elseif scenario == 2
        disp('Scenario 2: P12 with dynamic noise, increased by 2% (relative to attentiveness degree) per iteration')
        %Sample entries of P12 uniformly from [-0.5,0.5] and correct each 
        %row with a multiplication factor such that the sum of each the entries of each row sums up to 1
        P12_vec = rand(1,m2)-0.5;
        P12_norm = P12_vec / (P12_vec * ones(m2,1));
        P12 = ones(m1,1) * P12_norm;
 
        %Noise factor increased by 2% relative to attentiveness degree of
        %1 per iteration
        noise = 0.02 * i;
        %Sample random noise uniformly from [-noise, noise] for each entry
        P12noise = (rand(m1,m2)-0.5) * 2 * noise;
        
        %Add noise to P12; thus P12 is neither a 1-equal attentiveness
        %matrix or a matrix with constant columns anymore
        P12 = P12 + P12noise;
        
        %P21 is a random 1-equal attentiveness matrix with entries sampled
        %uniformly from [-0.5,0.5]
        P21 = rand(m2,m1)-0.5;
        P21_norm = inv(diag(P21 * ones(m1,1)));
        P21 = P21_norm * P21;
 
        %Define A
        A = [C(1,1)*P11 C(1,2)*P12; C(2,1)*P21 C(2,2)*P22]; 
        
        %SaveName = ['A_SensAA12Noise2PercentIteration' num2str(i)];
        %save(SaveName,'A');
        
        ExportDPI = 300;
        ExportName = 'SensAA12Noise2Percent300.png';
            
    elseif scenario == 3
        disp('Scenario 3: rows of groups with dynamics noise, increased by 0.1% (relative to attentiveness degree) per iteration')
        %Sample P12, a 1-equal attentiveness matrix with constant columns (as before)
        P12_vec = rand(1,m2)-0.5;
        P12_norm = P12_vec / (P12_vec * ones(m2,1));
        P12 = ones(m1,1) * P12_norm;
        
        %Multiply each row of P12 with a noise factor in sampled uniformly
        %from [1-0.001*i,1+0.001*i]
        dev = 0.001 * i;
        P12 = diag(1 + (rand(m1,1)-0.5) * dev) * P12;    
               
        %Construct P21 as a 1-equal attentiveness matrix with entries
        %randomly picked between [-0.5,0.5] before correction by a
        %multiplication factor
        P21 = rand(m2,m1)-0.5;
        P21_norm = inv(diag(P21 * ones(m1,1)));
        P21 = P21_norm * P21;
        %Add noise
        P21 = diag(1 + (rand(m2,1)-0.5) * dev) * P21;     
        
        %Construct P11 with the initial parameters
        P11 = [1.2 -0.2; 0.01 0.99];
        P11_norm = inv(diag(P11 * ones(m1,1)));
        P11 = P11_norm * P11;
        %Add noise
        P11 = diag(1 + (rand(m1,1)-0.5) * dev) * P11;
        
        %Construct P22 with the initial parameters
        P22 = [0.9 0.1; -0.2 1.2];
        P22_norm = inv(diag(P22 * ones(m2,1)));
        P22 = P22_norm * P22;
        
        %Add noise
        P22 = diag(1 + (rand(m2,1)-0.5) * dev) * P22;
        
        %Define A
        A = [C(1,1)*P11 C(1,2)*P12; C(2,1)*P21 C(2,2)*P22]; 
        
        %SaveName = ['A_SensAAttentNoise0_1PercentIteration' num2str(i)];
        %save(SaveName,'A');
        
        ExportDPI = 300;
        ExportName = 'SensAAttentNoise0_1Percent300.png';
        
    end

    %Calculation, sorting and saving of eigenvalues
    v_eigenvalues(:,i) = sort(abs(eig(A)));

end

figure
plot(1:1:steps,v_eigenvalues(:,:),'LineStyle','none','Marker','.','MarkerSize',8,'LineWidth',2)
xlabel('Iteration') 
ylabel('Eigenvalues of A') 
%give the figure the name f in preperation for export
f = gcf
%export with a resolution of ExportDPI and name ExportName
exportgraphics(f,ExportName,'Resolution',ExportDPI)