function output = v7Shocks(n,y,pReplace,pRecover,pAlly,allyEvolve,fitEvolve,varAll,varFit,cutoff,plots,start, distribution,both)
%v7Shocks(n,y,pReplace,pRecover,pAlly,allyEvolve,fitEvolve,varAll,varFit,cutoff,plots,start,distribution,both)
%Size of shock and number of agents
%affected by the shock are both picked from the chosen distribution
%'u' for uniform 'b' for beta 'e' for truncated exponential 'n' for normal
%if both=1 both shocksize and shock range are drawn from specified
%distribution - else only shock size is.
%n= # agents, y= # years,
%pRepace=probability of extinct agents being replaced,
%pRecover=prob of alive agents recovering health lost last period,
%pAlly=prob forming alliance,
%allyEvolve=whether or not alliances evolve,
%fitEvolve=whether or not fitness evolves,
%varAll=size of alliance variation,
%varFit=size of fitness variation,
%cutoff=minimum size of shock that is passed on,
%plot=whether or not to create otputs and plots,
%start=number of initial values to exclude from plots
if plots>0
tic
end
parameters.n=n;
parameters.y=y;
parameters.pReplace=pReplace;
parameters.pRecover=pRecover;
parameters.pAlly=pAlly;
parameters.allyEvolve=allyEvolve;
parameters.fitEvolve=fitEvolve;
parameters.varAll=varAll;
parameters.varFit=varFit;
parameters.cutoff=cutoff;
parameters.version='v7Shocks';
parameters.distribution=distribution;
parameters.both=both;
if plots>0
parameters
end
%Arrays
fitness=zeros(1,n);
Fitness=zeros(1,y);
bFitness=zeros(1,y);
shocksize=zeros(1,y);
shockcentre=0;
shocknumber=zeros(1,y);
Alive=zeros(1,y);
extinct=zeros(1,n);
extinctp=zeros(1,n);
Extinct=zeros(1,y);
Extinctp=zeros(1,y);
replaced=zeros(1,n);
Replaced=zeros(1,y);
recovered=zeros(1,n);
Recovered=zeros(1,y);
alliances=zeros(n);
Alliances=zeros(1,y);
bAlliances=zeros(1,y);
k=zeros(n,y);
bk=zeros(n,y);
hits=zeros(n);
lasthits=zeros(n);
dist=zeros(n);
position=sort([unifrnd(0,1,1,n)]);
for i=1:n
dist(i,:)=[min(abs(position-position(i)),abs(1-position+position(i)))];
end
%Year 1
initialfitness=unifrnd(0,1,1,n);
fitness=initialfitness;
lastfitness=initialfitness;
Alive(1)=n;
agents=1:n;
lastextinct=0;
shocked=0;
tmp=0;
%Start Loop
for t=2:y
extinct=zeros(1,n);
extinctp=zeros(1,n);
if plots>0 && round(t/100)==t/100
t
end
finish=t;
%Allow agents the chance to recover any health they just lost
if pRecover>0 && length(tmp)>0
recovered=zeros(1,n);
recovered(lastfitness~=fitness)=binornd(1,pRecover,1,length(lastfitness~=fitness));
Recovered(t)=sum(recovered);
fitness(recovered==1)=lastfitness(recovered==1);
end
%Now we give agents who went extinct last period the chance to be
%replaced
if pReplace>0 && sum(lastextinct)>0
replaced=zeros(1,n);
replaced(lastextinct)=binornd(1,pReplace,1,length(lastextinct));
fitness(replaced==1)=initialfitness(replaced==1);
agents=find(fitness>0);
end
if length(agents)==0
break
end
%If allyEvolve==1 subject every alliance to a variation - bound so
%fitness stays in (0,1)
if allyEvolve==1 && sum(sum(alliances>0))>0
diffm=zeros(n);
diffm(alliances>0)=unifrnd(-varAll,varAll,sum(sum(alliances>0)),1);
%find the total attempted change to each row and the scaling fact
total=sum(diffm');
scale=zeros(1,n);
scale(total>0)=1-fitness(total>0);
scale(total<0)=fitness(total<0);
%now find the new alliances value
alliances=diag(scale)*diffm+alliances;
fitness=sum((diag(scale)*diffm)')+fitness;
end
%If fitEvolve=1 all agents are subject to small changes in fitness
if fitEvolve==1
vari=zeros(1,n);
vari(agents)=unifrnd(-varFit,varFit,1,length(agents));
scale=zeros(1,n);
scale(vari>0)=1-fitness(vari>0);
scale(vari<0)=fitness(vari<0);
fitness=fitness+vari.*scale;
end
%Now allow new alliances to form between ALIVE agents
if pAlly>0
%find gaps where there are not currently alliances and applies a
%probabilty to choose which will ally
space=zeros(n);
space(agents,agents)=alliances(agents,agents)==0;
space=tril(binornd(1,pAlly,n),-1).*space;
if sum(sum(space))>0
[i,j]=find(space==1);
for s = 1:length(i)
v1=unifrnd(0,1);
v2=unifrnd(0,1);
alliances(i(s),j(s))=v1*(1-fitness(i(s)));
alliances(j(s),i(s))=v2*(1-fitness(j(s)));
fitness(i(s))=fitness(i(s))+alliances(i(s),j(s));
fitness(j(s))=fitness(j(s))+alliances(j(s),i(s));
end
end
end
%Count how many alive agents we have going in to the shock
Alive(t)=sum(fitness>0);
lastfitness=fitness;
bAlliances(t)=sum(sum(alliances>0))/2;
bFitness(t)=mean(fitness);
bk(:,t)=sum(alliances~=0)'; %find degree of each vertex
%Now we apply the shock!
%first pick out its size, centre, and the range over which the shock is
%applied. The matrix dist tells us the distance between agents.
%nb shock range must be in [0,0.5] since all locations are on
%a unit circumference circle.
shockrange(t)=unifrnd(0,1)/2;
if distribution=='e'
s=1;
while s==1
shocksize(t)=exprnd(0.2);
if shocksize(t)<=1
s=0;
end
end
if (both==1)
while s==1
shockrange(t)=exprnd(0.2);
if shockrange(t)<=1
s=0;
end
end
end
elseif distribution=='n'
s=1;
while s==1
shocksize(t)=normrnd(0.5,0.1);
if shocksize(t)<=1 && shocksize(t)>=0
s=0;
end
end
if (both==1)
while s==1
shockrange(t)=normrnd(0.5,0.1)/2;
if shockrange(t)<=0.5 && shockrange(t)>=0
s=0;
end
end
end
elseif distribution=='b'
shocksize(t)=betarnd(1,5);
if (both==1)
shockrange(t)=betarnd(1,5)/2;
end
else
shocksize(t)=unifrnd(0,1);
end
if length(agents)==1
shockcentre=agents;
else shockcentre=randsample(agents,1);
end
isshocked=dist(shockcentre,:)0);
for s=shocked
lasthits(s,s)=min(shocksize(t)/fitness(s), 1);
end
fitness(shocked)=fitness(shocked)-shocksize(t);
extinct(fitness<0)=1;
fitness(fitness<0)=0;
%now the shock propogates
while 1>0
agents=find(fitness>0);
if length(agents)==0
break
end
lasthits(lasthits0);
if fitness==0
scale=1;
else scale=1;
end
for j=from
to=agents(agents~=j);
hits(to,i)=lasthits(i,j)*scale*alliances(to,i);
end
end
%record hits as a % of fitness for next loop to use
tmp=fitness;
tmp(fitness==0)=1;
lasthits=min(hits*diag(1./tmp),1);
%now we have a matrix of the hits - so can remove dead alliances
lastalliances=alliances;
alliances=alliances-hits;
tmpm=zeros(n);
tmpm(agents,agents)=alliances(agents,agents);
alliances=tmpm;
%now apply the hits
fitness=fitness-sum(hits');
extinctp(fitness<0)=1;
fitness(fitness<0)=0;
end
alliances(fitness==0,:)=0;
alliances(:,fitness==0)=0;
%Almost at the end of the loop - record totals
Alliances(t)=sum(sum(alliances>0))/2;
k(:,t)=sum(alliances~=0)'; %find degree of each vertex
Extinct(t)=sum(extinct);
Extinctp(t)=sum(extinctp);
Fitness(t)=mean(fitness);
Recovered(t)=sum(recovered);
lastextinct=find((extinct+extinctp)>0);
%extinct=zeros(1,n);
%extinctp=zeros(1,n);
end
%Finished year loop - generate plots and outputs
%prepare outputs
output.parameters=parameters;
output.Alive=Alive;
output.Alliances=Alliances;
output.bAlliances=bAlliances;
output.Extinct=Extinct;
output.Extinctp=Extinctp;
output.Fitness=Fitness;
output.bFitness=bFitness;
output.Replaced=Replaced;
output.Recovered=Recovered;
output.shocksize=shocksize;
output.shocknumber=shocknumber;
output.finish=finish;
output.k=k;
output.bk=bk;
output.avek=mean(k);
output.bavek=mean(bk);
output.varavek=var(output.bavek);
output.bvaravek=var(output.bavek);
output.aveavek=mean(output.avek);
output.baveavek=mean(output.bavek);
output.varallk=var(reshape(k,1,n*y));
output.bvarallk=var(reshape(bk,1,n*y));
output.aveallk=mean(reshape(k,1,n*y));
output.baveallk=mean(reshape(bk,1,n*y));
if plots==1 | start>finish
%GRAPHS
figure
subplot(2,2,1)
plot(Alive(1:finish))
title('Alive')
subplot(2,2,2)
plot(Alliances(1:finish))
title('Alliances')
subplot(2,2,3)
scatter(1:finish, Extinct(1:finish),2)
title('Extinct')
subplot(2,2,4)
scatter(1:finish, Extinctp(1:finish),2)
title('Extinctp')
figure
subplot(2,2,1)
hist((Extinct(1:finish)+Extinctp(1:finish))./Alive(1:finish))
title('(Extinct(1:finish)+Extinctp(1:finish))./Alive(1:finish)')
subplot(2,2,2)
hist((Extinct(1:finish)+Extinctp(1:finish)))
title('(Extinct(1:finish)+Extinctp(1:finish))')
subplot(2,2,3)
scatter(1:finish, shocknumber(1:finish),2)
title('shocknumber')
subplot(2,2,4)
scatter(1:finish, shocksize(1:finish),2)
title('shocksize')
end
if plots==2 && start0)
toc
end