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(coop2, 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
		
find_best <- function(payoff, delta, phi,tol){
	nbest <- 1
	nequ <- 0
	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)
					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(welfa > mwelfx){
							mwelfx <- welfa
							mstratx <- c(xans[2], vcoop)
						}
					}
					if(mwelfx > -10000) nequ <- nequ + 1
					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))
	print(paste("nbest=", nbest, "nequ=",nequ))
	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
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)
					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), 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

ffac <- 1.0

#compute Delta cutoff
Delta <- function(pmat, xpi){
	xdelta <- pmat[1]- pmat[2]
	xdenom <- ((1-xpi)*(pmat[3] - pmat[1])) + (xpi * (pmat[4] - pmat[2]))
	return(xdelta/xdenom)
}

#compute beta cutoff
beta <- function(pmat, disc){
	U <- pmat[1]
	G <- pmat[3] - pmat[1]
	L <- pmat[4] - pmat[2]
	num <- ((1-disc)*G) - (disc * (pmat[1] - pmat[4]))
	den <- ((1-disc)*(G-L)) - (disc * (pmat[1] - pmat[4]))
	ans <- num / den
	if(num > 0) ans <- 0
	return(ans - (1/2))
}
	

#compute FK Delta cutoff
FK <- function(pmat, disc){
	xG <- pmat[3] - pmat[1]
	xUN <- pmat[1]- pmat[4]
	xL <- pmat[4] - pmat[2]
	dRD <- (xG + xL) / (xUN + xG + xL)
	return(disc - dRD)
}

#compute probabilities for different punishment strategies
LDD <- function(rho, xpi){
	x1 <- 1 - (2*xpi)
	x2 <- rho + (xpi * xpi)
	return( (x2/x1)/xpi )
}

dDD <- function(xpi){
	return(xpi * ( 1 - (2*xpi)))
}

pDD <- function(xpi){
	return(xpi * xpi)
}

dDC <- function(xpi){
		yy <- 1 - (2 * xpi)
		return(yy * yy)
}

pDC <- function(xpi){
	return(2 * xpi * (1 - xpi))
}

dall <-function(xpi){
	x1 <- 1 - xpi
	x2 <- 1 - (2 * xpi)
	return(x1 * x2)
}

pall <- function(xpi){
	return(xpi * ( 2 - xpi))
}

welfare <- function(rho, xpi, pmat, cdelta, pmixed){
	if(pmixed > 0) print("mixed upper bound")
	u1 <- (1 - xpi) * (1-xpi)* pmat[1]
	u2 <- xpi * (1-xpi) * pmat[3]
	u3 <- xpi * (1-xpi) * pmat[2]
	u4 <- xpi * xpi * pmat[4]
	u <- u1 + u2 + u3 + u4
	welf <- pmat[4] * (1 - xpi) * (1-xpi)
	welf <- welf + ((pmat[2]+pmat[3])* xpi * (1-xpi))
	welf <- welf + (pmat[1] * xpi *xpi)
	n <- welf
	u <- n + (ffac * (u-n))
	if(cdelta == -1) return(welf)
	G <- pmat[3] - pmat[1]
	L <- pmat[4] - pmat[2]
	xpii <- xpi + (1-ffac)
	g <- (1-2*xpii)*(((1-xpii) * G) + (xpii * L))
	ell <- (1-(2*xpi))*(((1-xpi)*L)+(xpi * G))
	if(pmixed > 0) g <- ((1-pmixed) * g )+ (pmixed * ell)
	if(LDD(rho, xpi) <= cdelta){
		print("constraint does not bind")
		dP <- dDD(xpi)
		xP <- pDD(xpi)
		welf <- xP/dP
		welf <- u - (welf * g)
		return(welf)
	}
	print("constraint does bind")
	p <- pDD(xpi)
	q <- p + dDD(xpi)
	p2 <- pDC(xpi)
	q2 <- p2 + dDC(xpi)
	rat1 <- cdelta * (q-p)/(rho + p)
	rat2 <- p2 * (rho + q)/(rho + p)
	denom <- q2 - rat2
	if(denom < 1/1000){
		print("totally not feasible")
		return(n)
	}
	P2 <- (1-rat1)/(q2 - rat2)
	P2 <- rho * g * P2
	print(paste("P2=", P2))

	vmn <- rho * (u - n) - (p2 * P2)
	vmn <- vmn /(rho + p)
	v <- n + vmn
	test <- v - P2 - n
	print(paste("test=", test))
	if(test < 0) {
		print("totally not feasible")
		return(n)
	}
	return(v)
}

arate <- function(xpi, pmat, welf){
	xCC <- pmat[1]
	DC <- pmat[3]
	CD <- pmat[2]
	DD <- pmat[4]
	u1 <- (1 - xpi) * (1-xpi)* xCC
	u2 <- xpi * (1-xpi) * DC
	u3 <- xpi * (1-xpi) * CD
	u4 <- xpi * xpi * DD
	u <- u1 + u2 + u3 + u4
	n <- DD * (1 - xpi) * (1-xpi)
	n <- n + ((DC+CD)* xpi * (1-xpi))
	n <- n + (xCC * xpi *xpi)
	lambda <- (welf - n)/(u - n)
	phi1 <- lambda * (1-xpi)*(1-xpi)
	phi2 <- (1-lambda)*xpi*xpi
	phi3 <- xpi * (1-xpi)
	phi <- phi1 + phi2 + phi3
	print(paste("theory accrate=", phi))
	
}
	
	
		
bestall <- function(xpi, disc, pmat){
			cdelta <- Delta(pmat, xpi)
			rho <- (1-disc)/disc
			nwelf <- welfare(rho, xpi, pmat, -1, 0)
			print(paste("nwelf=", nwelf))
			welf <- welfare(rho, xpi, pmat, cdelta, 0)
			arate(xpi, pmat, welf)
			mixrat <- (welf-pmat[4])/(pmat[1]-pmat[4])
			print(paste("theory ratio", (welf-pmat[4])/(pmat[1]-pmat[4])))
			pmixed <- 1
			if(mixrat > 0.5) {
				offf <- (pmat[2] + pmat[3])/2
				if(offf > 100) print(paste("offf=", offf))
				else pmixed <- sqrt(1-mixrat)
			}
			welf <- welfare(rho, xpi, pmat, cdelta, pmixed)
			mixrat <- (welf-pmat[4])/(pmat[1]-pmat[4])
			print(paste("mixed ratio", mixrat))
			
}


soptimize <- function(hh){
	zdata <- matis[which(matis$rep == hh),]
	ii <- get_treat(zdata, 1)
	experiment <- hbtreatments$experiment[ii]
	mat <- matrices[ii,]
	bestall(phi/2, deltas[ii], mat)
}
	

make_report <- function(hh, tol){
	ewelf <- get_welfare(hh)
	soptimize(hh)
	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
report_all()

