clear all;
%%% FINAL VERSION: 17 MARCH 2021 %%%
%%% THE FOLLOWING .MAT FILES ARE PART OF THIS CODE: (1)
%%% group_500_500_random1_per_st_op_A.mat and (2) group_500_500_random1_per_st_op_x0.mat

%C is chosen as an unequal attentiveness matrix, with (1.5 - 1)^2 < -4 *
%-0.3 * 0.5 such that C describes a system of two agents that exhibits
%asymptotic periodic behavior
C = [1.5 -0.3; 0.5 1];

%Choose two groups with sizes m1 and m2
m1 = 500;
m2 = 500;

%Entries are uniformly sampled between [-neg, 1-neg]
neg = 0.38;

%Entries of P11 are sampled uniformly from [-neg, 1-neg]
P11 = (rand(m1,m1)-neg)*2;

%Produce 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;

%Repeat the above steps for P12, P22 and P21
P12 = (rand(m1,m2)-neg)*2;
P12_norm = inv(diag(P12 * ones(m2,1)));
P12 = P12_norm * P12;

P22 = (rand(m2,m2)-neg)*2;
P22_norm = inv(diag(P22 * ones(m2,1)));
P22 = P22_norm * P22;

P21 = (rand(m2,m1)-neg)*2;
P21_norm = inv(diag(P21 * ones(m1,1)));
P21 = P21_norm * P21;

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

%For an exact copy of Figure 5 one should load the following A and x0 which
%were produced by following the steps above
%load('group_500_500_random1_per_st_op_A','A')
%load('group_500_500_random1_per_st_op_x0','x0')

%Store the number of agents as N
N = size(A,1);

%Create a new x0 with entries between [-0.5,0.5]
x0 = rand(N,1)-0.5;

%Number of iteration steps
steps = 50;

%Initiate a matrix of 0s to store the development of relative opinions in
vector = zeros(N,steps);

%Calculation of relative opinions with the chosen A and x0
for i = 1:1:steps
  if i == 1
    x = x0;
  else
    x = A * x;
  end
  %normalize such that the aggregate of the absolute entries of the vector equals
  %N
  x = x / (ones(1,N) * abs(x))*N;
  vector(:,i) = x;
end

%Plot relative opinions
figure
tiledlayout(3,3)
for i = 1:1:9
    nexttile
    grid on;
    xlim([-5, 5]);
    ylim([0,550]);
    title('Time')
    h1=histogram(vector(1:m1,i))
    h1.BinWidth=0.2;
    hold on
    h2=histogram(vector(m1+1:m1+m2,i))
    h2.BinWidth=0.2;
    title(['Time = ' num2str(i)])
    grid on;
    xlim([-5, 5]);
    ylim([0,550]);   
    hold off
end
%give the figure the name f in preperation for export
f = gcf
%export with a resolution of 300DPI
exportgraphics(f,'HistogramGroups300.png','Resolution',300)