# punt2.R

# library used
library(stringr)

# function used
mixmod.fun <- function(x0, mu0=10, eps=1e-8, maxiter=500) {
# this function takes a matrix with each row representing a different "sample"
# the first column is the observed frequency and the second column is the number of observations
# it outputs the posterior and prior distributions over the true probabilities associated
# with each sample.
  x0 <- as.matrix(x0)  
  p_hat <- x0[,1] 
  N_hat <- x0[,2]  
  N <- dim(x0)[1]
  if (length(mu0) > 1) {
    R <- length(mu0)
    mu <- mu0
  } else {
    R <- mu0 
    mu <- NULL
    for (r in 1:R) {
      mu <- c(mu, max((r-1)/R,0.0001))
    }
  }
  
  # Initial values
  z_hat <- matrix(1/R,N,R)
  z_init <- floor(p_hat*R) + 1
  z_init <- ifelse(z_init>R,R,z_init)
  for (i in 1:N) {
    z_hat[i,z_init[i]] <- 1
  }
  z_hat <- z_hat/rowSums(z_hat)
  iter <- 1
  finished <- FALSE
  lambda <- matrix(0, nrow = maxiter, ncol = R)
  lambda[1,] <- colMeans(z_hat)
  log_stuff <- lfactorial(N_hat) - lfactorial(p_hat*N_hat) - lfactorial(N_hat - p_hat*N_hat)
  Lf_hat_fixed <- matrix(NA,N,R)
  for (r in 1:R) {
    Lf_hat_fixed[,r] <- log_stuff + (p_hat*N_hat)*log(mu[r]) + (N_hat - p_hat*N_hat)*log(1 - mu[r])
  }
  
  while (!finished) {
    z_hat <- t(matrix(rep(lambda[iter,],N),ncol=N))
    f_hat <- matrix(NA,N,R)
    f_hat <- exp(Lf_hat_fixed)*z_hat
    Lf_hat <- Lf_hat_fixed + log(z_hat)
    for (i in 1:N) {
      if (sum(f_hat[i,],na.rm=TRUE)==0) {
        a0 <- order(-Lf_hat[i,])
        f_hat[i,a0[1]] <- 1
      }
    }
    z_hat <- f_hat/rowSums(f_hat,na.rm = TRUE)
    iter <- iter + 1
    lambda[iter,] <- colMeans(z_hat)
    finished <- iter >= maxiter
    maxchange <- sum((lambda[iter,] - lambda[iter-1,])^2)
    finished <- finished | (maxchange < eps)
    print(iter)
    print(maxchange)
  }
  return(list(data = x0, posteriors = z_hat, iter= iter, mu=mu, lambda = lambda[1:iter,], lambdahat = lambda[iter,]))
}


# Import data and create distance measures

x <- read.csv("../input/NFL-Punt-Analytics-Competition/play_information.csv", as.is = TRUE)
y <- read.csv("../input/NFL-Punt-Analytics-Competition/video_review.csv", as.is = TRUE)
x$concussion <- NA
x$concussion <- ifelse(paste(x$GameKey,x$PlayID,sep = "") %in% unique(paste(y$GameKey,y$PlayID,sep = "")),1,0)
x$dist <- as.numeric(sub(" ","",sub(".[A-z]*","",x$YardLine)))
x$side <- sub("[ ].*","",x$YardLine)
x$dist <- ifelse(x$side==x$Poss_Team,100-x$dist,x$dist)

# This step takes a while
# An alternative is to use a pre-made data:
x <- read.csv("../input/nfl-punt-data-with-kick-distance/punt.csv", as.is = TRUE)
# x$punt_dist <- NA
# for (i in 1:length(strsplit(x$PlayDescription," "))) {
#   for (j in 1:length(strsplit(x$PlayDescription," ")[[i]])) {
#     if (strsplit(x$PlayDescription," ")[[i]][j]=="punts") {
#       x[i,]$punt_dist <- as.numeric(strsplit(x$PlayDescription," ")[[i]][j+1])
#     }
#   }
#  print(i)
# }


# Effect of Rule
# probability of concussion in opp half
sum(ifelse(x[x$concussion==1,]$dist < 50,1,0), na.rm = TRUE)/sum(ifelse(x$dist < 50,1,0), na.rm = TRUE)
# probability of concussion in own half
sum(ifelse(x[x$concussion==1,]$dist > 49,1,0), na.rm = TRUE)/sum(ifelse(x$dist > 49,1,0), na.rm = TRUE)
# probability of concussion in own half short punt
sum(ifelse(x[x$concussion==1,]$dist > 49 & x[x$concussion==1,]$punt_dist < 40,1,0), na.rm = TRUE)/sum(ifelse(x$dist > 49 & x$punt_dist < 40,1,0), na.rm = TRUE)
# probability of concussion in own half long putn
sum(ifelse(x[x$concussion==1,]$dist > 49 & x[x$concussion==1,]$punt_dist > 39,1,0), na.rm = TRUE)/sum(ifelse(x$dist > 49 & x$punt_dist > 39,1,0), na.rm = TRUE)




# Create Concussion Probabilities in 10 x 10 yards
# of distance to goal (dist) and punt distance (punt_dist)
# it calculates the number of punts and the number of concussions for each to go distance x punt distance.

eps <- 0.00001
prob <- matrix(NA,8*7,6)
min_d <- 30 - eps
i <- 1
for (d in c(4:10)) {
  max_d <- d*10 
  min_p <- 0 - eps
  for (p in c(1:8)) {
    max_p <- p*10 
    prob[i,1] <- sum(ifelse(x$dist > min_d & x$dist < max_d & x$punt_dist > min_p & x$punt_dist < max_p,1,0), na.rm = TRUE)
    prob[i,2] <- sum(ifelse(x[x$concussion==1,]$dist > min_d & x[x$concussion==1,]$dist < max_d & x[x$concussion==1,]$punt_dist > min_p & x[x$concussion==1,]$punt_dist < max_p,1,0), na.rm = TRUE)
    prob[i,3:6] <- c(min_d+eps,max_d,min_p+eps,max_p)
    min_p <- max_p - eps
    i <- i + 1
    print(i)
  }
  min_d <- max_d - eps
}


p1 <- prob[,2]/prob[,1]  # calculates the raw probabilities
N1 <- prob[,1]  # the number of punts
V1 <- prob[,3:6]  # mid-point yardage.
p2 <- p1[is.nan(p1)==0]  # removes the "areas" with no punts.
N2 <- N1[is.nan(p1)==0]
V2 <- V1[is.nan(p1)==0,]


# this function takes the raw frequencies and uses them to calculate the posterior probabilities of 
# a concussion an each 10 x 10 area.
a <- mixmod.fun(cbind(p2,N2), mu0 = 10000, eps = 1e-12, maxiter = 100000)

a_post_mean <- a$posteriors%*%as.vector(c(1:10000))/10000
mean_a_post_mean <- mean(a_post_mean)

res <- cbind(V2,a_post_mean)  # this collects the mean posterior for each 10 x 10 (to go distance x punt distance)



# 20 x 20  
# same as above but for a bigger area.

eps <- 0.00001
prob3 <- matrix(NA,4*4,6)
min_d <- 30 - eps
i <- 1
for (d in c(1:4)) {
  max_d <- 30 + d*20 
  min_p <- 0 - eps
  for (p in c(1:4)) {
    max_p <- p*20 
    prob3[i,1] <- sum(ifelse(x$dist > min_d & x$dist < max_d & x$punt_dist > min_p & x$punt_dist < max_p,1,0), na.rm = TRUE)
    prob3[i,2] <- sum(ifelse(x[x$concussion==1,]$dist > min_d & x[x$concussion==1,]$dist < max_d & x[x$concussion==1,]$punt_dist > min_p & x[x$concussion==1,]$punt_dist < max_p,1,0), na.rm = TRUE)
    prob3[i,3:6] <- c(min_d+eps,max_d,min_p+eps,max_p)
    min_p <- max_p - eps
    i <- i + 1
    print(i)
  }
  min_d <- max_d - eps
}


p13 <- prob3[,2]/prob3[,1]
N13 <- prob3[,1]
V13 <- prob3[,3:6]
p23 <- p13[is.nan(p13)==0]
N23 <- N13[is.nan(p13)==0]
V23 <- V13[is.nan(p13)==0,]

d <- mixmod.fun(cbind(p23,N23), mu0 = 10000, eps = 1e-11, maxiter = 100000)

d_post_mean <- d$posteriors%*%as.vector(c(1:10000))/10000

res_d <- cbind(V23,d_post_mean)