% GNU Public Licence Copyright (c) Colin Torney
% Comments and questions to colin.j.torney@gmail.com
% This code is provided freely, however when using this code you are asked to cite our related paper:
% Colin J. Torney, Tommaso Lorenzi, Iain D. Couzin, Simon A. Levin (2015) Social information use and the evolution of unresponsiveness in collective systems
clear all
close all
clc
wG=0.2;
NA=200;
k=4;
NT=1000;
tau=100; % (period of the external signal)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% external signal (E(t))
t=linspace(1,NT,NT);
E=(2/pi)*asin(sin((2*pi/tau)*(t)));
averaging = 100;
% start all agents with the same wS
for wS = 0.5:0.01:0.6
SuccessRate=0;
for av = 1:averaging
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G=zeros(NT,NA);
u=zeros(NT,NA);
AVSign=zeros(1,NT);
ap=zeros(NT,NA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initial data
G(1,:)=0;
u(1,:)=(2*randi([0 1],1,NA))-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% create all the G for all times by using an OU process around E, then
% calculate a from the integration of the Gaussian distribution around Gi
for i=1:NA
a = exp(-wG);
b = E(1:NT-1)*(1-a);
c = sqrt((1-a*a)/2/wG);
G(:,i) = [0 filter(1,[1 -a], b+c*randn(1,NT-1), 0)]';
ap(:,i) = (erf((sqrt(wG))*(1+G(:,i))) - erf((sqrt(wG))*G(:,i)))./(erf((sqrt(wG))*(1-G(:,i))) + erf((sqrt(wG))*G(:,i)));
end
for p=2:NT
%% rewire the network
pos_edges = zeros(NA-1,NA);
a1 = zeros(NA,NA);
% exclude self from list of possible neighbours
for i=1:NA
offset = 0;
for j=1:NA-1
if (j==i) offset = 1; end
pos_edges(j,i) = j+offset;
end
end
% randomly create k distinct connections
for i=1:NA
a1(i,pos_edges(randperm(NA-1,k),i))=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%evaluate \Deltan
np=zeros(1,NA);
nn=zeros(1,NA);
for i=1:NA
%% calculate the social force for one side or the other
for j=1:NA
if a1(i,j)>0
%\Deltan(u_i)
if u(p-1,j)>0
np(i)=np(i)+1;
else
nn(i)=nn(i)+1;
end
end
end
% calculate the difference in opinions
deltan=np(i) - nn(i);
% weight according to the social confidence level
PR = ap(p,i)*(((1 - wS)/wS)^deltan);
%Evaluate u_i(t)
if abs(PR-1) < 1e-6
u(p,i)=(2*randi([0 1],1))-1;
else
u(p,i) = sign(1-PR);
end
end
AVSign(p) = 1/NA*sum(u(p,:)==sign(E(p)));
end
SuccessRate=SuccessRate + 1/NT*sum(AVSign);
end
fprintf('%f %f\n',wS, SuccessRate/averaging);
end