rm(list=ls())
dyn.load("/n/nobackup2/ellison_lab/srecord/src/glm_amcmc.so")
source("/n/nobackup2/ellison_lab/srecord/src/glm_amcmc.R")
source("/n/nobackup2/ellison_lab/srecord/src/mkMatUtil.R")
library(coda)

##Read data
treebp0 <- read.csv("/n/nobackup2/ellison_lab/srecord/spGLM/spenv_alltimesteps.csv", header=TRUE)

# subsample data maintaining 10% of the data for a hold-out validation
treebp0 <- subset(treebp0, tenths<=8)

tsuga.pa <- treebp0[,5]
fagus.pa <- as.matrix(treebp0[,4])
weights <- as.matrix(treebp0[,3])
x <- as.matrix(treebp0[,c(6,11)])

##Collect samples
t_fit <- glm((tsuga.pa/weights)~x, weights=weights, family="binomial")
summary(t_fit)

f_fit <- glm((fagus.pa/weights)~x, weights=weights, family="binomial")
summary(f_fit)

t_beta.starting <- coefficients(t_fit)
t_beta.tuning <- t(chol(vcov(t_fit)))

f_beta.starting <- coefficients(f_fit)
f_beta.tuning <- t(chol(vcov(f_fit)))

n.batch <- 500
batch.length <- 200
n.samples <- n.batch*batch.length

# Run model for tsuga
t_nonspatialc1 <- glm.amcmc(tsuga.pa~x, family="binomial", weights=weights, 
                 starting=list("beta"=t_beta.starting),
                 tuning=list("beta"=t_beta.tuning),
                 priors=list("beta.flat"),
                 amcmc=list("n.batch"=n.batch,"batch.length"=batch.length,"accept.rate"=0.43),
                 n.samples=n.samples, sub.samples=c(100,n.samples,10),
                 verbose=TRUE, n.report=10)

#pdf(file="/n/nobackup2/ellison_lab/srecord/spGLM/tsuga_nonspatial.pdf")
#plot(t_nonspatial$p.samples)
#dev.off()

print(summary(t_nonspatialc1$p.samples))
print(t_nonspatialc1$DIC)
save(t_nonspatialc1, file="/n/nobackup2/ellison_lab/srecord/spGLM/nonspGLM_tsugac1.Rdata")

# Run model for fagus
f_nonspatialc1 <- glm.amcmc(fagus.pa~x, family="binomial", weights=weights, 
                 starting=list("beta"=f_beta.starting),
                 tuning=list("beta"=f_beta.tuning),
                 priors=list("beta.flat"),
                 amcmc=list("n.batch"=n.batch,"batch.length"=batch.length,"accept.rate"=0.43),
                 n.samples=n.samples, sub.samples=c(100,n.samples,10),
                 verbose=TRUE, n.report=10)

#pdf(file="/n/nobackup2/ellison_lab/srecord/spGLM/fagus_nonspatialc1.pdf")
#plot(f_nonspatial$p.samples)
#dev.off()

print(summary(f_nonspatialc1$p.samples))
print(f_nonspatialc1$DIC)

save(f_nonspatialc1, file="/n/nobackup2/ellison_lab/srecord/spGLM/nonspGLM_fagusc1.Rdata")