% Code to generate data simulating the voting experiments using parameters
% provided by the fitting methods.

clear all

tic

% define groups sizes
N_vals = [12,13,16,18,19,20,22,25,28,31,33];

% Pick one specifically
% N_vals = [31];

for Nit = 1:length(N_vals)

    N = N_vals(Nit);
    
    % read the fitted parameters
    params = readmatrix(['./fitted_rates/rates_',num2str(N),'_1.txt'])
    
    arl = params(1);
    brl = params(2);
    crl = params(3);
    ars = params(4);
    brs = params(5);
    crs = params(6);
    asr = params(7);
    bsr = params(8);
    csr = params(9);
    
    % time discretisation (1 round)
    TD = 1;
    
    % duration of the simulation (number of rounds)
    T = 5*60*60*10;

    % arrays to record votes
    votes = zeros(1,N);
    votes(1) = 1;
    rec_votes = zeros(N,T);
    rec_votes(1:N,1) = votes;
    
    for t=1:T-1
            % find state of the system
            nr = sum(votes==0); % votes for X
            nl = sum(votes==1); % votes for Y
            
            % computes transition rate matrix of the model
            Q = zeros([3 3 length(nr)]);
            Q(1,2) = arl + brl.*nl/N + crl.*(N-nr-nl)/N;
            Q(1,3) = ars + brs.*(N-nr-nl)/N + crs.*nl/N;
            Q(1,1) = -Q(1,2)-Q(1,3);
            Q(2,1) = arl + brl.*nr/N + crl.*(N-nr-nl)/N;
            Q(2,3) = ars + brs.*(N-nr-nl)/N + crs.*nr/N;
            Q(2,2) = -Q(2,1)-Q(2,3);
            Q(3,1) = asr + bsr.*nr/N + csr.*nl/N;
            Q(3,2) = asr + bsr.*nl/N + csr.*nr/N;
            Q(3,3) = -Q(3,1)-Q(3,2);

            % computes transition matrix (probabilities of transitioning
            % between states)
            eQ = expm(Q*TD);

            % loop through each agent
            for j=1:N
                % finds the vote of the agent in the previous round
                pv = votes(j);
                % finds the probability that the agent transitions to each
                % state
                vr = eQ(pv+1,1); % prob of transition to X
                vl = eQ(pv+1,2); % prob of transition to Y
                vabs = eQ(pv+1,3); % prob of abstaining
                coin = rand; % flip a coin
                % now determines which transition happened
                if coin < vr
                    votes(j)=0;
                elseif coin < vr+vl
                    votes(j)=1;
                else
                    votes(j)=2;
                end
            end
            % Record the votes
            rec_votes(:,t+1)=votes;
    end
    % gets the aggregate number of votes for each option
    vts(1,:) = sum(rec_votes==0); % X
    vts(2,:) = sum(rec_votes==1); % Y
    Vs = [vts(1,:)',vts(2,:)']; % rewrites as a matrix

     % saving the data
    fname = sprintf("./simulations/simdata_%d_1.txt",N);
    % writematrix(Vs,fname); % uncomment to save the data
end

% plot time series
T=Vs;
figure(1)
plot(1:length(Vs(:,1)), Vs)

% plot hsitogram
figure(2)
histogram(Vs(:,1)-Vs(:,2))

toc