mono <- read.csv("monotone.22.csv")

nseed <- 40538

inite <- 0
alpha <- -0.268
beta <- 1.291
lambda <- 0.182
pcc <- 0.995
pcd <- 0.355
pdd <- 0.012

g <- mono$DC - 1
l <- -mono$CD
Drd <-(g+l)/(1+g+l)
Drd <- mono$delta - Drd

get0 <- function(Drd, e){
	p0 <- alpha + (beta * Drd) + e
	p0 <- 1 + exp(-p0)
	p0 <- 1/p0
	return(p0)
}

updatee <- function(lage, coop, V){
	coeff <- -1
	if(coop == 1) coeff <- 1
	enew <- lambda * coeff * V
	enew <- enew + lage
	return(enew)
}

get_payoff <- function(coop1, coop2, CD, DC){
#coop is zero for defect, one for cooperate
	pay <- c(1, CD, DC, 0)
	ind <- (2 * (1-coop1)) + (2 - coop2)
	return(pay[ind])
}

getC <- function(pr){
	coop <- 0
	x <- runif(1,0.0,1.0)
	if(x <= pr) coop <- 1
	return(coop)
}

rcont <- function(delta){
	cont <- 0
	x <- runif(1,0.0,1.0)
	if(x <= delta) cont <- 1
	return(cont)
}

nexto <- function(c1, c2){
	y <- pcd
	if((c1 * c2) == 1){
		#both cooperated
		return(pcc)
	}
	if((c1 + c2) == 0){
		#both defected
		return(pdd)
	}
	return(y)
}	
	

supersim <- function(inite1, inite2, Drd, CD, DC, delta){
	V1 <- 0
	V2 <- 0
	coop1 <- 0
	coop2 <- 0
	round <- 0
	continue <- 1
	p1 <- get0(Drd, inite1)
	p2 <- get0(Drd, inite2)
	while(continue == 1){
		c1 <- getC(p1)
		c2 <- getC(p2)
		round <- round + 1
		if(round == 1){
			coop1 <- c1
			coop2 <- c2
		}
		v1 <- get_payoff(c1, c2, CD, DC)
		v2 <- get_payoff(c2, c1, CD, DC)
		V1 <- V1 + v1
		V2 <- V2 + v2
		p1 <- nexto(c1,c2)
		p2 <- nexto(c2, c1)
		continue <- rcont(delta)
	}	
	welf <- (V1 + V2)/2
	welf <- welf/round
	e1 <- updatee(inite1, coop1, V1)
	e2 <- updatee(inite2, coop2, V2)
	return(list(welf, e1, e2))
}

set.seed(nseed)

ngroups <- 7
npop <- 2 * ngroups

fetch_groups <- function(){
	groups <- matrix(0, ngroups, 2)
	glist <- list()
	for(ii in 1:ngroups){
		for(jj in 1:npop){
			if(!(jj %in% glist)){
				groups[ii,1] <- jj
				glist <- append(glist, jj)
				continue <- 1
				while(continue == 1){
					rn = floor(runif(1, min=ii + 1, max=npop + 1))
					if(!(rn %in% glist)){
						groups[ii,2] <- rn
						glist <- append(glist, rn)
						continue <- 0
					}
				}
				break
			}
		}
	}
	return(groups)
}
		
sim <- function(ii){
	statex <- rep(0,npop)
	Drd <- Drd[ii]
	CD <- mono$CD[ii]
	DC <- mono$DC[ii]
	delta <- mono$delta[ii]
	matches <- mono$matches[ii]
	welf <- 0
	nrec <- 0
	for(ii in 1:matches){
		#fetch a ngroup x 2 matrix of matches
		groups <- fetch_groups()
		for(jj in 1:ngroups){
			ans <- supersim(statex[groups[jj,1]], statex[groups[jj,1]], Drd, CD, DC, delta)
			if(ii >= 10){
				welf <- welf + ans[[1]]
				nrec <- nrec + 1
			}
			statex[groups[jj,1]] <- ans[[2]]
			statex[groups[jj,2]] <- ans[[3]]
		}
	}
	welf <- welf/nrec
	return(welf)
}

iter <- 1000

multi_sim <- function(ii){
	welf <- 0
	for(jj in 1:iter){
		welf <- welf + sim(ii)
	}
	return(welf/iter)
}

for(ii in 1:nrow(mono)){
	print(paste(mono$experiment[ii], mono$delta[ii], multi_sim(ii)))
}				



