require("lattice")
require("MCMCpack")

# extract credible intervals for the dist.fullributions
dist.full <- read.csv("dist.full.csv",header=T,stringsAsFactors=F)

if (0) { # ignore for now

# output file
ofile <- file("dist.full.mean_ci.csv","w")
cat("season,location,pstart,",sep="",file=ofile);
cat(colnames(dist.full[4:NCOL(dist.full)]),sep=",,",file=ofile)
cat("\n",file=ofile)

for (i in 1:NROW(dist.full)) {
	cat(dist.full[i,1],",",dist.full[i,2],",",dist.full[i,3],sep="",file=ofile)
	for (j in 4:NCOL(dist.full)) {
		alpha <- 1+as.numeric(dist.full[i,j])
		beta <- 0
		for (k in 4:NCOL(dist.full)) {
			if (j != k) { beta <- beta+1+as.numeric(dist.full[i,k]) }
		}

		qs <- qbeta(c(0.5,0.025,0.975), alpha, beta)
		if (qs[1] < 0.01) {
			cat(",",sprintf("%.4f%%",qs[1]*100),",\"(",sprintf("%.4f%%",qs[2]*100),", ",sprintf("%.4f%%",qs[3]*100),")\"",sep="",file=ofile)
		} else {
			cat(",",sprintf("%.1f%%",qs[1]*100),",\"(",sprintf("%.1f%%",qs[2]*100),", ",sprintf("%.1f%%",qs[3]*100),")\"",sep="",file=ofile)
		}
	}
	cat("\n",file=ofile)
}

close(ofile)

} # ignore for now

# create images

# this function draws figure based on inputs
"draw_fig" <- function(away,home,away_min,away_max,home_min,home_max,ylimits,title,xaxis,yaxis,xnames,fname) {
	# setup panel
	panel.q <- function (x, y, ...) {
	  ltyp <- "solid"
	  lwidth <- 1

		lpoints(x=1:NROW(home)+0.15,y=as.numeric(home),pch=19,col="blue")
		for (i in 1:NROW(away_min)) {
			llines(x=c(i,i)-0.15,y=as.numeric(c(away_min[i],away_max[i])),col="red",lty=ltyp,lwd=lwidth)
			llines(x=c(i,i)+0.15,y=as.numeric(c(home_min[i],home_max[i])),col="blue",lty=ltyp,lwd=lwidth)
		}

	  panel.xyplot(x, y, ...)
	}

	# create Y labels
	ylab <- c()
	for (i in ylimits[1]:ylimits[2]) {
		if (i %% ylimits[3] == 0) {
			if (yaxis == "Points") {
				ylab <- c(ylab,paste(i,"",sep=""))
			} else {
				ylab <- c(ylab,paste(i,"%",sep=""))
			}
		} else {
			ylab <- c(ylab,"")
		}
	}

	if (yaxis == "Points") {
		ylims <- range(ylimits[1],ylimits[2])
	} else {
		ylims <- range(30,65)
	}

	fig <- xyplot(as.numeric(away) ~ 1:NROW(away)-0.15, type="p", col="red", pch=19, panel=panel.q,
		main=title, xlab=xaxis, ylab=yaxis,
		#ylim=range(c( as.numeric(c(away_min,home_min))-0.5, as.numeric(c(away_max,home_max))+0.5 )),
		#ylim=range(30,65),
		ylim=ylims,
		scales=list(x=list(at=seq(1,NROW(xnames),1),labels=xnames,rot=45),y=list(at=seq(ylimits[1],ylimits[2],ylimits[4]),labels=ylab))
	)

	# save image as PDF
	trellis.device(dev=pdf, file=fname); #, width="5", height="3"); # default width=height=7

	print(fig);

	# save the file
	graphics.off();
}


# These images use basic rules: miss+SF is not an attempt; make+SF is an attempt
# 2FG%
items <- c(1,2,4,6,8,9,10,12,14,17,18,19,20)
items <- c(items,items+20)

away <- c()
home <- c()
away_min <- c()
away_max <- c()
home_min <- c()
home_max <- c()
xnames <- c()

for (i in items) {
	observations <- sum(dist.full[c(i,i+40,i+80),4:NCOL(dist.full)])
	fg2_miss <- sum(dist.full[c(i,i+40,i+80),c(5,8)])
	fg2_make <- sum(dist.full[c(i,i+40,i+80),c(6,7)])

	probs <- rdirichlet(10000,c(1+fg2_miss,1+fg2_make,1+observations-fg2_miss-fg2_make))
	cprobs <- probs[,2]/(1-probs[,3]) # probability of 2FGM given a 2FGA
	qs <- sprintf("%.1f", quantile(cprobs,c(0.5,0.025,0.975))*100)

	if (dist.full[i,2] == "away") {
		away <- c(away,qs[1])
		away_min <- c(away_min,qs[2])
		away_max <- c(away_max,qs[3])
		xnames <- c(xnames,dist.full[i,3])
	} else {
		home <- c(home,qs[1])
		home_min <- c(home_min,qs[2])
		home_max <- c(home_max,qs[3])
	}

#	if ( (as.numeric(qs[3])-as.numeric(qs[2])) <= 5) {
#cat(i,"\n")
#cat(i,": ",dist.full[i,1],"/",dist.full[i,2],"/",dist.full[i,3],": ",qs[1],"% (",qs[2],"%, ",qs[3],"%)\n",sep="")
#	}

}

draw_fig(away,home,away_min,away_max,home_min,home_max,c(1,100,5,1),"06-07 to 08-09 2FG% Given Start of Play","Play Start","2FG%",xnames,"2fg_pct.pdf");

# 3FG%
items <- c(1,2,4,6,8,9,10,12,14,17,18,19,20)
items <- c(items,items+20)

away <- c()
home <- c()
away_min <- c()
away_max <- c()
home_min <- c()
home_max <- c()
xnames <- c()

for (i in items) {
	observations <- sum(dist.full[c(i,i+40,i+80),4:NCOL(dist.full)])
	fg3_miss <- sum(dist.full[c(i,i+40,i+80),c(10,13)])
	fg3_make <- sum(dist.full[c(i,i+40,i+80),c(11,12)])

	probs <- rdirichlet(10000,c(1+fg3_miss,1+fg3_make,1+observations-fg3_miss-fg3_make))
	cprobs <- probs[,2]/(1-probs[,3]) # probability of 3FGM given a 3FGA
	qs <- sprintf("%.1f", quantile(cprobs,c(0.5,0.025,0.975))*100)

	if (dist.full[i,2] == "away") {
		away <- c(away,qs[1])
		away_min <- c(away_min,qs[2])
		away_max <- c(away_max,qs[3])
		xnames <- c(xnames,dist.full[i,3])
	} else {
		home <- c(home,qs[1])
		home_min <- c(home_min,qs[2])
		home_max <- c(home_max,qs[3])
	}

#	if ( (as.numeric(qs[3])-as.numeric(qs[2])) <= 5) {
#cat(i,"\n")
#cat(i,": ",dist.full[i,1],"/",dist.full[i,2],"/",dist.full[i,3],": ",qs[1],"% (",qs[2],"%, ",qs[3],"%)\n",sep="")
#	}

}

draw_fig(away,home,away_min,away_max,home_min,home_max,c(1,100,5,1),"06-07 to 08-09 3FG% Given Start of Play","Play Start","3FG%",xnames,"3fg_pct.pdf");

# eFG%
items <- c(1,2,4,6,8,9,10,12,14,17,18,19,20)
items <- c(items,items+20)

away <- c()
home <- c()
away_min <- c()
away_max <- c()
home_min <- c()
home_max <- c()
xnames <- c()

for (i in items) {
	observations <- sum(dist.full[c(i,i+40,i+80),4:NCOL(dist.full)])
	fg2_miss <- sum(dist.full[c(i,i+40,i+80),c(5,8)])
	fg2_make <- sum(dist.full[c(i,i+40,i+80),c(6,7)])
	fg3_miss <- sum(dist.full[c(i,i+40,i+80),c(10,13)])
	fg3_make <- sum(dist.full[c(i,i+40,i+80),c(11,12)])

	probs <- rdirichlet(10000,c(1+fg2_miss,1+fg2_make,1+fg3_miss,1+fg3_make,1+observations-fg2_miss-fg2_make-fg3_miss-fg3_make))
	# eFG% x 2 = points per shot attempt iff eFG% = points per shot attempt / 2
	cprobs_fg2 <- probs[,2]/(1-probs[,5]) # probability of 2FGM given a shot attempt
	cprobs_fg3 <- probs[,4]/(1-probs[,5]) # probability of 3FGM given a shot attempt
	cprobs <- cprobs_fg2 + (cprobs_fg3 * 3)/2
	qs <- sprintf("%.1f", quantile(cprobs,c(0.5,0.025,0.975))*100)

	if (dist.full[i,2] == "away") {
		away <- c(away,qs[1])
		away_min <- c(away_min,qs[2])
		away_max <- c(away_max,qs[3])
		xnames <- c(xnames,dist.full[i,3])
	} else {
		home <- c(home,qs[1])
		home_min <- c(home_min,qs[2])
		home_max <- c(home_max,qs[3])
	}

#	if ( (as.numeric(qs[3])-as.numeric(qs[2])) <= 5) {
#cat(i,"\n")
#cat(i,": ",dist.full[i,1],"/",dist.full[i,2],"/",dist.full[i,3],": ",qs[1],"% (",qs[2],"%, ",qs[3],"%)\n",sep="")
#	}

}

draw_fig(away,home,away_min,away_max,home_min,home_max,c(1,100,5,1),"06-07 to 08-09 eFG% Given Start of Play","Play Start","eFG%",xnames,"efg_pct.pdf");

# Points per Play: Assuming Pr(Make FT) = 0.75 (uses all 2FG and 3FG shots)
# This tells us how many points to expect given a shot attempt, as it takes
# free throw shots into account
items <- c(1,2,4,6,8,9,10,12,14,17,18,19,20)
items <- c(items,items+20)

away <- c()
home <- c()
away_min <- c()
away_max <- c()
home_min <- c()
home_max <- c()
xnames <- c()

for (i in items) {
	observations <- sum(dist.full[c(i,i+40,i+80),4:NCOL(dist.full)])
	fg2_miss <- sum(dist.full[c(i,i+40,i+80),c(5,8)])
	fg2_make <- sum(dist.full[c(i,i+40,i+80),c(6)])          # excludes Make+SF
	fg3_miss <- sum(dist.full[c(i,i+40,i+80),c(10,13)])
	fg3_make <- sum(dist.full[c(i,i+40,i+80),c(11)])         # excludes Make+SF

	fg2_miss_sf <- sum(dist.full[c(i,i+40,i+80),c(9)])
	fg2_make_sf <- sum(dist.full[c(i,i+40,i+80),c(7)])
	fg3_miss_sf <- sum(dist.full[c(i,i+40,i+80),c(14)])
	fg3_make_sf <- sum(dist.full[c(i,i+40,i+80),c(12)])

	probs <- rdirichlet(10000,c(1+fg2_miss,1+fg2_make,1+fg3_miss,1+fg3_make,1+fg2_miss_sf,1+fg2_make_sf,1+fg3_miss_sf,1+fg3_make_sf,1+observations-fg2_miss-fg2_make-fg3_miss-fg3_make-fg2_miss_sf-fg2_make_sf-fg3_miss_sf-fg3_make_sf))
	cprobs_fg2 <- probs[,2]/(1-probs[,9]) # probability of 2FGM given a shot event
	cprobs_fg3 <- probs[,4]/(1-probs[,9]) # probability of 3FGM given a shot event
	cprobs_fg2_miss_sf <- probs[,5]/(1-probs[,9]) # probability of 2 FT
	cprobs_fg3_miss_sf <- probs[,7]/(1-probs[,9]) # probability of 3 FT
	cprobs_fg2_make_sf <- probs[,6]/(1-probs[,9]) # probability of 2 points + 1 FT
	cprobs_fg3_make_sf <- probs[,8]/(1-probs[,9]) # probability of 3 points + 1 FT
	cprobs <- cprobs_fg2*2 + cprobs_fg3 * 3 +
		cprobs_fg2_miss_sf*(2*0.75) + cprobs_fg3_miss_sf*(3*0.75) +
		cprobs_fg2_make_sf*(2.75) + cprobs_fg3_make_sf*(3.75)
	qs <- sprintf("%.2f", quantile(cprobs,c(0.5,0.025,0.975))*100)

	if (dist.full[i,2] == "away") {
		away <- c(away,qs[1])
		away_min <- c(away_min,qs[2])
		away_max <- c(away_max,qs[3])
		xnames <- c(xnames,dist.full[i,3])
	} else {
		home <- c(home,qs[1])
		home_min <- c(home_min,qs[2])
		home_max <- c(home_max,qs[3])
	}

#	if ( (as.numeric(qs[3])-as.numeric(qs[2])) <= 5) {
#cat(i,"\n")
#cat(i,": ",dist.full[i,1],"/",dist.full[i,2],"/",dist.full[i,3],": ",qs[1]," (",qs[2],", ",qs[3],")\n",sep="")
#	}

}

draw_fig(away,home,away_min,away_max,home_min,home_max,c(90,140,10,1),"06-07 to 08-09 Points per 100 Shot Events Given Start of Play","Play Start","Points",xnames,"pts_per_shot_event.pdf");