options(stringsAsFactors=FALSE)

#first observation
sobs <- 10

cpayoff <- c(100, 50, 50, 0)

get_joint <- function(coop1, coop2, phi){
	joint <- rep(0, 4)
	coop1 <- ((1-(phi/2)) * coop1) + ((phi/2) * (1 - coop1))
	coop2 <- ((1-(phi/2)) * coop2) + ((phi/2) * (1 - coop2))
	joint[1] <- coop1 * coop2
	joint[2] <- coop1 * (1-coop2)
	joint[3] <- (1-coop1) * coop2
	joint[4] <- (1-coop1) * (1-coop2)
	return(joint)
}

	
get_first <- function(payoff, coop1, coop2, phi){
	joint <- get_joint(coop1, coop2, phi)
	first <- 0
	for(ii in 1:4){
		first <- first + (joint[ii] * payoff[ii])
	}
	return(first)
}

get_values <- function(payoffs, vcoop1, vcoop2, delta, phi){
	A <- matrix(0, 4, 4)
	for(y in 1:4){
		A[y,] <- get_joint(vcoop1[y], vcoop2[y], phi)
	}	
	A <- diag(c(1,1,1,1)) - (delta * A)
	first <- rep(0, 4)
	for (y in 1:4){
		first[y] <- get_first(payoffs, vcoop1[y], vcoop2[y], phi)
	}
	first <- (1 - delta) * first
	return(solve(A, first))
}

get_first <- function(payoff, coop1, coop2, phi){
	joint <- get_joint(coop1, coop2, phi)
	first <- 0
	for(ii in 1:4){
		first <- first + (joint[ii] * payoff[ii])
	}
	return(first)
}

get_fpmat <- function(payoff, V, delta, phi){
	gmat <- rep(0, 4)
	for(ii in 1:4){
		coop1x <- 1
		coop2x <- 1
		if(ii == 2){
			coop1x <- 1
			coop2x <- 0
		}
		if(ii == 3){
			coop1x <- 0
			coop2x <- 1
		}
		if(ii == 4){
			coop1x <- 0
			coop2x <- 0
		}
		nf <- get_first(payoff, coop1x, coop2x, phi)
		jointx <- get_joint(coop1x, coop2x, phi)
		gmat[ii] <- (1 - delta) * nf
		for(yp in 1:4) gmat[ii] <- gmat[ii] + (delta * jointx[yp] * V[yp])
	}
	return(gmat)
}

get_mixed <- function(c,d){
	if((c >= 0) &&  (d >= 0)) return(0)
	if((c <= 0) && (d <= 0)) return(1)
	coop <- d/(d-c)
	return(coop)
}

get_max <- function(c,d, tol){
	if(d >= 0){
		md <- d + tol
		mc <- c + tol
		if(mc > 0) mc <- 0
		coop <- get_mixed(mc, md)
		return(coop)
	}
	md <- d - tol
	mc <- c - tol
	if(mc < 0) mc <- 0
	coop <- get_mixed(mc, md)
	return(coop)
}

get_min <- function(c,d, tol){
	if(d >= 0){
		md <- d - tol
		mc <- c - tol
		if(md < 0) md <- 0
		coop <- get_mixed(mc, md)
		return(coop)
	}
	md <- d + tol
	mc <- c + tol
	if(md > 0) md <- 0
	coop <- get_mixed(mc, md)
	return(coop)
}

get_eq_pay <- function(gmat, coop){
	eq_pay1 <- (coop^2) * gmat[1]
	eq_pay2 <-coop * (1-coop) * (gmat[2] + gmat[3])
	eq_pay4 <- ((1-coop)^2) * gmat[4]
	return(eq_pay1 + eq_pay2 + eq_pay4)
} 

get_best_first <- function(gmat, b1, b2, b3, b4){
	sderiv <- 2 * (gmat[1] - gmat[2] - gmat[3] + gmat[4])
	small <- b1
	if(b2 < b1) small <- b2
	large <- b4
	if(b3 > b4) large <- b3
	xbests <- small
	xbest <- get_eq_pay(gmat, small)
	xbest2 <- get_eq_pay(gmat, large)
	if(xbest2 > xbest) {
		xbest <- xbest2
		xbests <- large
	}
	if(sderiv >= 0){
		return(c(xbest, xbests))
	}

	p <- gmat[2] + gmat[3] - (2 * gmat[4])
	p <- -p
	p <- p/sderiv
	if((p >= b2) && (p <= b3)) {
		best <- get_eq_pay(gmat, p)
		if(xbest > best) return(c(xbest, xbests))
		return(c(best,p))
	}
	bv <- c(b1, b2, b3, b4)
	tt <- rep(0, 4)
	for(ii in 1:4) tt[ii] <- abs(bv[ii] - p)
	bbt <- 1
	bb <- tt[1]
	for(ii in 2:4){
		if(tt[ii] < bb) {
			bbt < ii
			bb <- tt[ii]
		}
	}
	bests <- bv[bbt]
	best <- get_eq_pay(gmat, bests)
	if(xbest > best) return(c(xbest, xbests))
	return(c(best, bests))
}	

get_eq <- function(gmat, payoff, V, delta, phi, tol){
	c <- gmat[3] - gmat[1]
	d <- gmat[4] - gmat[2]
	cc <- 0
	dd <- 0
	if(c > tol) cc <- 1
	if(c < -tol) cc <- -1
	if(d > tol) dd <- 1
	if(d < -tol) dd <- -1
	#return domC, domD, coord, anti, all, low_mixed, high_mixed
	if((cc == -1) && (dd == -1)) return(c(1, 0, 0, 0, 0, 1, 1))
	if((cc == 1) && (dd == 1)) return(c(0, 1, 0, 0, 0, 0, 0))
	if((cc == 0) && (dd == 0)) return(c(0, 0, 0, 0, 1, 0, 1))
	minc <- get_min(c, d, tol)
	maxc <- get_max(c, d, tol)
	if(cc == 1){
		if(dd == -1) return(c(0, 0, 0, 1, 0, minc, maxc))
		return(c(0, 1, 0, 1, 0, minc, maxc))
	}
	if(cc == -1){
		if(dd == 1) return(c(0, 0, 1, 0, 0, minc, maxc))
		return(c(1, 0, 1, 0, 0, minc, maxc))
	}
	if(dd == 1) return(c(0, 1, 1, 0, 0, minc, maxc))
	return(c(1, 0, 0, 1, 0, minc, maxc))	
}

gel <- function(neq, base, ans){
	if(neq == 1) amat <- ans
	else amat <- cbind(base, ans)
	return(amat)
}

check_eq <- function(coop1, coop2, payoff, V, delta, phi, tol){
	gmat <- get_fpmat(payoff, V, delta, phi)
	ans <- get_eq(gmat, payoff, V, delta, phi, tol)
	#domC, domD, coord, anti, all, low_mixed, high_mixed
	if(ans[5] == 1) return(get_best_first(gmat, 0, 0, 1, 1))
	is_eq <- 0
	bestv <- c(0,0)
	for(ii in 1:4) {
		if(ans[ii] == 1){
			if(ii == 1){
				if((coop1 == 1) && (coop2 == 1)) {
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, 1, 1, 1, 1))
				}
			}
			if(ii == 2){
				if((coop1 == 0) && (coop2 == 0)){
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, 0, 0, 0, 0))
				}
			}
			if(ii > 2){
				lcoop <- min(coop1, coop2)
				hcoop <- max(coop1, coop2)
				if((lcoop >= ans[6]) && (hcoop <= ans[7])) {
					if(coop1 == coop2){
						is_eq <- is_eq + 1
						bestv <- gel(is_eq, bestv, get_best_first(gmat, ans[6], ans[6], ans[7], ans[7]))
					}
				}
			}
			if(ii == 3){
				if((coop1 == 1) && (coop2 == 1)) {
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, 1, 1, 1, 1))
				}
				if((coop1 == 0) && (coop2 == 0)) {
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, 0, 0, 0, 0))
				}
			}
			if(ii == 4){
				#in the anti-symmetric case the only symmetric first period equilibrium is the mixed one, checked above
				#this checks for the pure strategy asymmetric ones
				if((coop1 == 1) && (coop2 == 0)) {
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, ans[6], ans[6], ans[7], ans[7]))
				}
				if((coop1 == 0) && (coop2 == 1)) {
					is_eq <- is_eq + 1
					bestv <- gel(is_eq, bestv, get_best_first(gmat, ans[6], ans[6], ans[7], ans[7]))
				}
			}
		}
	}

	if(is_eq == 0) return(c(-10000, 0))
	themax <- -10000
	themaxv <- c(-10000, 0)
	for(ii in 1:is_eq){
		if(is_eq == 1) gbf <- bestv
		else gbf <- bestv[,ii]
		if(gbf[1] > themax) {
			themax <- gbf[1]
			themaxv <- gbf
		}
	}
	return(themaxv)
}

#final run	
inc <- 80
inc0 <- inc
inc0 <- 1

#base run; commented out for final run
#inc <- 20
#inc0 <- 20

#coarse run
#inc <- 10

#for bru_kam only
#inc <- 160
		
find_best <- function(payoff, delta, phi,tol){
	nbest <- 1
	nequ <- 0
	vequ <- c()
	bot <- 0
	vcoop <- rep(0, 4)
	mwelfare <- -100000
	mstrat <- rep(0,5)
	for(ii in 0:inc0){
		vcoop[1] <- ii/inc0
		for(jj in bot:inc){
			vcoop[2] <- jj/inc
			for(kk in bot:inc){
				vcoop[3] <- kk/inc
				for(ll in bot:inc){
					vcoop[4] <- ll/inc
					vcoop2 <- vcoop
					vcoop2[2] <- vcoop[3]
					vcoop2[3] <- vcoop[2]
					V <- get_values(payoff, vcoop, vcoop2, delta, phi)
					mwelfx <- -1000
					mstratx <- rep(0,5)
					lwelf <- -10000
					hwelf <- -10000
					for(hhh in 1:4) {
						if(mwelfx < -10000) break
						xans <- check_eq(vcoop[hhh], vcoop2[hhh], payoff, V, delta, phi, tol)
						welfa <- xans[1]
						if(welfa < -1000) {
							mwelfx <- -100000
							break
						}
						if(hhh == 1) lwelf <- welfa
						if(hhh == 4) hwelf <- welfa
						if(welfa > mwelfx){
							mwelfx <- welfa
							mstratx <- c(xans[2], vcoop)
						}
					}
					if(mwelfx > -10000) {
						nequ <- nequ + 1
						vequ <- append(vequ, lwelf)
						if(hwelf != lwelf){
							nequ <- nequ + 1
							vequ <- append(vequ, hwelf)
						}
					}
					if(mwelfx == mwelfare) nbest <- nbest + 1
					if(mwelfx > mwelfare){
						nbest <- 1
						mwelfare <- mwelfx
						mstrat <- mstratx
					}
				}
			}
		}
	}

	vcoop1 <- rep(0,4)
	for(ii in 1:4) vcoop1[ii] <- mstrat[ii + 1]
	vcoop2 <- vcoop1
	vcoop2[2] <- vcoop1[3]
	vcoop2[3] <- vcoop1[2]
	V <- get_values(payoff, vcoop1, vcoop2, delta, phi)
	cV <- get_values(cpayoff, vcoop1, vcoop2, delta, phi)
	crate <-cV[4] 
	if(mstrat[1] == 1) crate <- cV[1]
	gmat <- get_fpmat(payoff, V, delta, phi)
	print(gmat)
	print(get_eq (gmat, payoff, V, delta, phi, tol))
	nclose <- 0
	for(llmm in 1:NROW(vequ)){
		vequ[llmm] <- relative_welfare(vequ[llmm], payoff)
		if((relative_welfare(mwelfare, payoff) - vequ[llmm]) < 0.02) nclose <- nclose + 1
	}
	vequ <- sort(vequ)
	print(paste("nbest=", nbest, "nequ=",nequ, "nclose=", nclose))
	print(vequ)
	print(paste("theory coop=", relative_welfare(crate, cpayoff)))
	return(list(mwelfare, mstrat))
}



hbtreatments <- read.csv("hbtreatments.csv")
hbmetadata <- read.csv("hbmetadata.csv")
hbmatrices <- read.csv("hbmatrices.csv")
hbdata <- read.csv("hbdata.csv")

fktreatments <- read.csv("fktreatments.csv")
fkmetadata <- read.csv("fkmeta.csv")
fkmatrices <- read.csv("fkmatrices.csv")
fkdata <- read.csv("fkpdata.csv")

hbtreatments <- rbind(hbtreatments, fktreatments)
hbmetadata <- rbind(hbmetadata, fkmetadata)
hbmatrices <- rbind(hbmatrices, fkmatrices)
hbdata <- rbind(hbdata, fkdata)

hbmetadata <- hbmetadata[which(hbmetadata$game == "prisoner's dilemma"),]
exps <- unique(hbmetadata$experiment)
hbtreatments <- hbtreatments[which(hbtreatments$experiment %in% exps),]
hbtreatments <- hbtreatments[which(hbtreatments$supergame_ending < 1),]
hbdata <- hbdata[which(hbdata$match_period >= sobs),]
#hbdata <- hbdata[which(hbdata$match_period <= 30),]

exps <- unique(hbtreatments$experiment)
hbmetadata <- hbmetadata[which(hbmetadata$experiment %in% exps),]

row_num <- rep(0,nrow(hbtreatments))
hbtreatments <- cbind(hbtreatments, row_num)
for(ii in 1:nrow(hbtreatments)) hbtreatments$row_num[ii] <- ii

matis <- read.csv("mats.csv")

get_matrix <- function(experiment, payoff_treatment){
	matrix_row <- hbmatrices[which(hbmatrices$experiment == experiment),]
	matrix_row <- matrix_row[which(matrix_row$payoff_treatment == payoff_treatment),]
	matrix <- rep(0,4)
	for(ii in 1:4) {
		matrix[ii] <- matrix_row[1, paste("X", ii, sep = "")]
	}
	bot <- matrix[4]
	for(ii in 1:4) matrix[ii] <- matrix[ii] - bot
	defl <- matrix[1]
	matrix <- matrix/defl
	matrix <- 100 * matrix
	return(matrix)
}

get_matrices <- function(){
	matrices <- matrix(0, nrow(hbtreatments), 4)
	for(ii in 1:nrow(hbtreatments)){
		matrices[ii,] <- get_matrix(hbtreatments$experiment[ii], hbtreatments$payoff_treatment[ii])
	}
	return(matrices)
}

get_deltas <- function(){
	deltas <- rep(0, nrow(hbtreatments))
	for(ii in 1:nrow(hbtreatments)){
		deltas[ii] <- hbtreatments$supergame_ending[ii]
	}
	return(deltas)
}

get_treat <- function(udata, ii){
	vdata <- hbtreatments[which(hbtreatments$treatment == udata$treatment[ii]),]
	vdata <- vdata[which(vdata$experiment == udata$experiment[ii]),]
	return(vdata$row_num[1])
}

matrices <- get_matrices()
deltas <- get_deltas()

get_welfare <- function(hh){

zdata <- matis[which(matis$rep == hh),]
welf <- 0
coopr <- 0
obs <- 0

ccond <- c(0,0,0,0)
cctot <- c(0,0,0,0)
ccc1 <- 0
cct1 <- 0

nplay <- 0
nmatches <- 0
nvalid <- 0
tses <- 0
matrix <- rep(0,4)
wses <- list()

firstwelf <- 0
nfirst <- 0

for(kk in 1:nrow(zdata)){
	ii <- get_treat(zdata, kk)
	#each match appears twice in new data, so each observation only counts 1/2
	pfact <- 1
	if(zdata$fr[kk] == 1) pfact <- 1/2
	experiment <- hbtreatments$experiment[ii]
	treatment <- hbtreatments$treatment[ii]
	matrix <- matrices[ii,]
	udata <- hbdata[which(hbdata$experiment == experiment),]
	udata <- udata[which(udata$treatment == treatment),]
	part <- unique(udata$player)
	npart <- NROW(part)
	nplay <- nplay + (pfact * npart)
	sessions <- unique(udata$session)
	tses <- tses + NROW(sessions)

	for(jj in 1:NROW(sessions)){
		sesswelf <- 0
		sessno <- 0
		vdata <- udata[which(udata$session == sessions[jj]),]
		matches <- unique(vdata$match_period)
		if(NROW(matches) < 5) return(-2001)
		nmatches <- nmatches + NROW(matches) + 9
		nvalid <- nvalid + 1
		for(kk in 1:NROW(matches)){
			wdata <- vdata[which(vdata$match_period == matches[kk]),]
			groups <- unique(wdata$match_group)
			for(ll in 1:NROW(groups)){
				xdata <- wdata[which(wdata$match_group == groups[ll]),]
				rounds <- unique(xdata$supergame_period)
				for(mm in 1:NROW(rounds)){
					ydata <- xdata[which(xdata$supergame_period == rounds[mm]),]
					a1 <- ydata$a1[1]
					a2 <- ydata$a1[2]
					aa <- (2*(a2 - 1)) + a1
					bb <- (2*(a1 - 1)) + a2
					thiswelf <- (matrix[aa] + matrix[bb])/2
					welf <- welf + (pfact * thiswelf)
					if(kk == 1){
						if(mm == 1){
							firstwelf <- firstwelf + (pfact *thiswelf)
							nfirst <- nfirst + pfact
						}
					}
					sesswelf <- sesswelf + (pfact * thiswelf)
					coopr <- coopr + (pfact * (2 - a1 + 2 - a2)/2)
	
					ccround <- rounds[mm]
					if(ccround > 1){
						ccdata <- xdata[which(xdata$supergame_period == (ccround - 1)),]
						ca1 <- ccdata$a1[1]
						ca2 <- ccdata$a1[2]
						ccind <- (2 * (ca1 - 1)) + ca2
						if(a1 == 1) ccond[ccind] <- ccond[ccind] + (pfact/2)
						occind <- ccind
						if(occind == 2) {
							ccind <- 3
						}							
						if(occind == 3) {
							ccind <- 2
						}
						if(a2 == 1) ccond[ccind] <- ccond[ccind] + (pfact/2)
						cctot[occind] <- cctot[occind] + (pfact/2)
						cctot[ccind] <- cctot[ccind] + (pfact/2)
					}
					else{
						if(a1 == 1) ccc1 <- ccc1 + (pfact/2)
						if(a2 == 1) ccc1 <- ccc1 + (pfact/2)
						cct1 <- cct1 + pfact
					}
						

					obs <- obs + pfact
					sessno <- sessno + pfact
				}
			}
		}
		wses <- append(wses, sesswelf/(sessno * 100))
	}
}
print(paste(zdata$experiment[1], zdata$delta[1], matrix[2], matrix[3],  "ewelf", relative_welfare(welf/obs,matrix), "coop", coopr/obs, "nplay", nplay, "nmatches", nmatches/nvalid, "nplay/session", nplay/tses, "nsessions", tses), quote = FALSE)
print(ccc1/cct1)
print(ccond/cctot, quote = FALSE)
print(cctot/sum(cctot), quote = FALSE)
print(cctot/(sum(cctot)+cct1), quote = FALSE)

ws <- unlist(wses)
med <- median(ws)
top <- ws[which(ws >= med)]
bot <- ws[which(ws < med)]
if(NROW(ws) > 1) print(paste("min/max session welf", mean(bot), mean(top)), quote = FALSE)
else print("one session", quote = FALSE)

return(list(nplay, tses))
}

relative_welfare <- function(welfare, mat){
	rel <- welfare - mat[4]
	diff <- mat[1] - mat[4]
	return(rel/diff)
}

phi <- 1/3

optimize <- function(hh, tol){
	zdata <- matis[which(matis$rep == hh),]
	ii <- get_treat(zdata, 1)
	experiment <- hbtreatments$experiment[ii]
	mat <- matrices[ii,]

	#if(experiment != "bru_kam") return()
	#if(experiment != "dre_et_al_08") return()
	#if(experiment != "she_tar_sai"){
	#	if(experiment != "kag_sch") return()
	#}
	#if(hbtreatments$supergame_ending[ii] != 0.75) return()
	#if(hh != 10) return()


	best <- find_best(mat, deltas[ii], phi,tol)
	twelf <- best[[1]]
	mstrat <- best[[2]]
	retwelf <- relative_welfare(twelf, mat)
	print(paste("theory welf", retwelf, "strategy below"), quote = FALSE)
	print(mstrat, quote = FALSE)
	mmstrat <- (1-phi) * mstrat
	mmstrat <- mmstrat + (phi/2) * c(1,1,1,1,1)
	print(mmstrat, quote = FALSE)
	print("", quote = FALSE)
}


make_report <- function(hh, tol){
	ewelf <- get_welfare(hh)
	
	optimize(hh,tol)
	return(ewelf)
}

report_all <- function(){
	rest <- unique(matis$rep)
	nplay <- 0
	tses <- 0
	for(ii in 1:NROW(rest)) {
		ewelf <- make_report(ii,tol)
		nplay <- nplay + ewelf[[1]]
		tses <- tses + ewelf[[2]]
	}
	print(paste("av per sess", nplay/tses))
}

tol <- 0.2
#double tolerance
#tol <- 0.4
#halve tolerance
#tol <- 0.1
#quadruple tolerance
#tol <- 0.8

#phi <- 0.337

report_all()


for(ii in 0:10){
	#phi <- ii/10
	phi <- 0.2 + (ii/50)
	print(phi)
	report_all()
}

phi <- 0.332
report_all()

phi <- 0.333
report_all()

phi <- 0.338
report_all()

phi <- 0.339 
report_all()

phi <- 0.35
report_all()


phi = 1/3
omatrices <- matrices
matrices[25,3] <- omatrices[25,3] + 2
report_all()
matrices <- omatrices


