# Evolutionary Games in Stochastic Environments

library(MASS)

###############################
# Evolutionary games function #
###############################

# this function assumes there are a finite number of random variables determining
# the payoff matrix and that these random variables are all log-normally distributed

# INPUT
# n - number of strategies
# logvar.mean the log-means of the variables
# logvar.cov the covariance matrix of the log variables
# payoff a function of the variables that produces the payoff matrix
# d fraction that die
# T length of run
# x0 initial condition

# OUTPUT: x all frequencies over all time; rows correspend to different times

EG=function(n,logvar.mean,logvar.cov,payoff,T,x0,d){
	x=matrix(0,T,n)
	x[1,]=x0
	for(t in 1:(T-1)){
		var=exp(mvrnorm(1,mu=logvar.mean,Sigma=logvar.cov))
		A=payoff(var)
		x[t+1,]=(1-d)*x[t,]+d*(diag(x[t,])%*%A%*%x[t,])/sum(x[t,]*(A%*%x[t,]))
	}
	return(x)
}


###############
# SAMPLE USES #
###############

################################
# Lottery model with n species #
################################

# INPUTS:
# n number of competing species
# fitness.min minimum log fitness of the species
# fitness.max maximum log fitness of the species
# sigma log-variance in fitness (same for all)
# rho correlation in log-fitnesses between species
# d fraction dying  
# T length of simulation

Lottery=function(n=2,fitness.min=1,fitness.max=1.1,sigma=0,rho=1,d=0.9,T=1000){
	logvar.mean=seq(fitness.min,fitness.max,length=n) 	
	logvar.cov=(1-rho)*sigma*diag(n)+rho*sigma # the covariance matrix
	payoff=function(var)var%o%rep(1,n) # the payoff function
	x0=rep(1/n,n) # initially equal frequencies of all types

	x=EG(n=n,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d) #calling the command

	par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1)
	matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data
	par(mar=c(4,0,1,1))
	hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="")	
}


##################
# Hawk dove game #
##################

# INPUTS:
# V.mean mean payoff for a victory
# sigma log-variance in payoff for a victory
# C cost of a fight
# Base base payoff for everyone
# d fraction dying  
# T length of simulation



HawkDove=function(V.mean=2,C=2,Base=3,sigma=0,d=0.1,T=5000){
	logvar.cov=sigma
	logvar.mean=log(V.mean)-sigma/2 
	payoff=function(x)matrix(c((x-C)/2,x,
							0,x/2),2,2,byrow=TRUE)+Base
	x0=c(1/2,1/2)
	x=EG(n=2,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d)
	par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1)
	matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data
	par(mar=c(4,0,1,1))
	hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="")}
	
	
######################
# Rock Paper Scissor #
######################

# INPUTS:
# B.mean mean payoff for a win
# B.sigma log-variance for a win
# rho correlation between strategy payoffs for a win
# C cost for a loss (needs to less than base payoff)
# a.mean mean base payoff for everyone
# a.sigma log-variance base payoff for everyone
# rho correlation between log-wins and log-base
# d fraction dying
# T length of simulation



RPS=function(B.mean=0.5,B.sigma=0,rho=0,C=1,a.mean=2,a.sigma=0,d=0.1,T=5000){
  logvar.cov=matrix(0,6,6)
  logvar.cov[1:3,1:3]=diag(c(B.sigma,B.sigma,B.sigma))
  logvar.cov[4:6,4:6]=diag(c(a.sigma,a.sigma,a.sigma)) 
  logvar.cov[1:3,4:6]=diag(rep(rho*sqrt(B.sigma*a.sigma),3))
  logvar.cov[4:6,1:3]=diag(rep(rho*sqrt(B.sigma*a.sigma),3))
  logvar.mean=c(rep(log(B.mean)-B.sigma/2,3),rep(log(a.mean)-a.sigma/2,3))  
  payoff=function(x)matrix(c(x[4],x[4]-C,x[4]+x[1],x[5]+x[2],x[5],x[5]-C,x[6]-C,x[6]+x[3],x[6]),3,3,byrow=TRUE)
  x0=c(0.4,0.3,0.3)
  x=EG(n=3,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d)
  par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1)
  matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data
  par(mar=c(4,0,1,1))
  hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="")
}