################################################################################
# Program: SECR_Wojan_MultiSession_v7_1.R
# Directory: C:\MabryWojanSECR\
#
# Programmer: Shannon Knapp
# Date: 8 Dec 2013 (original)
# Updated: 20 Dec 2013
# Updated: 13 June 2014, 
# Updated: 14 June 2014 (renamed to "v5")
# Updated: 24 June 2014 - version 6, modified from version 5
# Updated: 27 June 2014 - version 7
# Updated: 6 July 2014 - Version 7.1
#
# Description:  Have multiple trapping sessions of data over 3 trapping grids
#
# Version 2b (builds off version 1, not version 2):
#	instead of only using the more northern 2 trapping grids, 
#	I will use all 3 grids, but have two disjoint sampling spaces 
#	S2 = S-North and S1 = S-South
#	
#	Also bringing in the dispersal points as in Version 3
#
# Version 3b (builds off version 2b)
#	adding on comparison of change in landscape-wide density
#
# Version 5 added the regression of distance against natal density
#
# Version 6 is being set-up to run on the cloud and to run BOTH sets of
#	dispersers simultaneously (with separate statistics for each)
#	But for ONLY one local/neighborhood radius
#
# Data Formatting:  
#	See C:\MabryWojanSECR\
#	Knapp20131018.msg & Wojan20131105.msg
#
# Code based on Gardner, Royle, & Wegan 2009. Ecology 90(4): 1106-1115
#	Assumes indiv caught in at most 1 trap per night, multiple trap nights
#
###############################################################################

# rm(list=ls())

###############################################################################
### Load libraries and set working directory
###############################################################################

library(ggplot2)
library(R2jags)

### Working directory
 dir <- "C:/MabryWojanSECRcloud"
# dir <- "C:/MabryWojanSECR"
setwd(dir)

### Model name
model.name <- "SECR_Wojan_MultiSession_v7_1"

date.today <- format(Sys.Date(), format="%Y%m%d")

### Workspace Filename
file.ws <- paste(model.name,"_",date.today,".RData", sep="")
# load(file.ws)

### Output (text) file name
outfile <- paste(model.name,"_",date.today,"_out.txt", sep="")

###############################################################################
### Ist his run a Trial/test run (MCMC) or Actual?  
###	1 = trial run
###	2 = actual run
###############################################################################

RunType <- 2

###############################################################################
### Which Disperser Data Set?    
###	1 = the smaller/telemetry set
###	2 = the larger/trapping set
###	3 = both at the same time
###############################################################################

DisperserSet <- 3		### Must be 1 or 2 or 3


###############################################################################
### Import data
###############################################################################

### Trap locations
trapcoords <- read.table(paste(dir,"/MultiSession20131112traps.txt",sep=""),header=TRUE) 

### H-matrix
Hmatrix <- read.table(paste(dir,"/MultiSession20131112captures.txt",sep=""),header=TRUE) 

### Session information
SessionInfo <- read.table(paste(dir,"/MultiSession20131112sessions.txt",sep=""),header=TRUE) 

### Dispersal (Natal and Settle locations)
if(DisperserSet == 1){
	### If Set 1
		DispWide <- read.csv(paste(dir,"/NatalSettleLocations20131118.csv",sep=""),header=TRUE)
		DispLong <- read.csv(paste(dir,"/NatalSettleLocations20131118long.csv",sep=""),header=TRUE)
		DispWide.original <- read.csv(paste(dir,"/NatalSettleLocations20131118.csv",sep=""),header=TRUE)
}else if(DisperserSet == 2){
	### If Set 2
		DispWide <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103.csv",sep=""),header=TRUE)
		DispLong <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103long.csv",sep=""),header=TRUE)
		DispWide.original <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103.csv",sep=""),header=TRUE)
}else if(DisperserSet == 3){
		DispWide.1 <- read.csv(paste(dir,"/NatalSettleLocations20131118.csv",sep=""),header=TRUE)
		DispLong.1 <- read.csv(paste(dir,"/NatalSettleLocations20131118long.csv",sep=""),header=TRUE)
		DispWide.original.1 <- read.csv(paste(dir,"/NatalSettleLocations20131118.csv",sep=""),header=TRUE)

		DispWide.2 <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103.csv",sep=""),header=TRUE)
		DispLong.2 <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103long.csv",sep=""),header=TRUE)
		DispWide.original.2 <- read.csv(paste(dir,"/NatalSettleLocationsSet2_20140103.csv",sep=""),header=TRUE)
		
		DispWide <- rbind(DispWide.1, DispWide.2)
		DispLong <- rbind(DispLong.1, DispLong.2)
		DispWide.original <-rbind(DispWide.original.1, DispWide.original.2)
}


###############################################################################################################
### Calculate disperal Distance
###############################################################################################################
DispWide$distance <- ((DispWide$natalX - DispWide$settleX)^2 + (DispWide$natalY - DispWide$settleY)^2)^0.5
DispWide.original$distance <- ((DispWide$natalX - DispWide$settleX)^2 + (DispWide$natalY - DispWide$settleY)^2)^0.5

 

###############################################################################################################
### Plot traps
###############################################################################################################
#ggplot(trapcoords, aes(x=UTMx, y=UTMy)) + 
#	geom_point(alpha=0.5, colour="red", shape=2, size=1) +
#	scale_shape(solid = FALSE) +
#	theme(panel.background=element_rect(fill="transparent", colour="NA")) +
#	coord_fixed()+
#	geom_line(data=DispLong, aes(x=UTMx,y=UTMy, group=ID), colour="blue")
#

###############################################################################
### Label traps as North=2 (traps 25-72) vs. South=1 (traps 1-24)
###############################################################################
#hist(trapcoords$UTMy)
trapcoords$group[trapcoords$UTMy >= 4262500] <- 2
trapcoords$group[trapcoords$UTMy < 4262500] <- 1

Hmatrix$group[Hmatrix$trap <= 24] <- 1
Hmatrix$group[Hmatrix$trap > 24] <- 2

###############################################################################################################
### Modify Plot Coordinates and set buffer
###############################################################################################################

BufferSize <- 75

trapcoords$Xnew <- trapcoords$UTMx - min(trapcoords$UTMx) + BufferSize
trapcoords$Ynew <- trapcoords$UTMy - min(trapcoords$UTMy) + BufferSize

S1x.min <- min(trapcoords$Xnew[trapcoords$group==1]) - BufferSize
S1x.max <- max(trapcoords$Xnew[trapcoords$group==1]) + BufferSize
S1y.min <- min(trapcoords$Ynew[trapcoords$group==1]) - BufferSize
S1y.max <- max(trapcoords$Ynew[trapcoords$group==1]) + BufferSize

S2x.min <- min(trapcoords$Xnew[trapcoords$group==2]) - BufferSize
S2x.max <- max(trapcoords$Xnew[trapcoords$group==2]) + BufferSize
S2y.min <- min(trapcoords$Ynew[trapcoords$group==2]) - BufferSize
S2y.max <- max(trapcoords$Ynew[trapcoords$group==2]) + BufferSize

Area.S1.SqMeters <- (S1x.max-S1x.min)*(S1y.max-S1y.min)
Area.S2.SqMeters <- (S2x.max-S2x.min)*(S2y.max-S2y.min)
Area.S.SqMeters <- Area.S1.SqMeters + Area.S2.SqMeters
Area.S.Hectares <- Area.S.SqMeters/10000
Area.S1.Hectares <- Area.S1.SqMeters/10000
Area.S2.Hectares <- Area.S2.SqMeters/10000
Area.S.Hectares
Area.S1.Hectares
Area.S2.Hectares
RegionAreaHectares <- c(Area.S1.Hectares, Area.S2.Hectares)


###############################################################################################################
### Plot traps
###############################################################################################################
#ggplot(trapcoords, aes(x=Xnew, y=Ynew)) + 
#	geom_point(alpha=0.5, colour="red", shape=2, size=1) +
#	scale_shape(solid = FALSE) +
#	theme(panel.background=element_rect(fill="transparent", colour="NA")) +
#	coord_fixed() +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.min,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.min,y=S1y.max,xend=S1x.max,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.max,y=S1y.max,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.min,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.min,y=S2y.max,xend=S2x.max,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.max,y=S2y.max,xend=S2x.max,yend=S2y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.max,yend=S2y.min)) 
#

###############################################################################################################
### Modify coordinates for dispersers and make sure they are IN S (natal and settle locations) 
###############################################################################################################

### DispBuffer -- in meters, the radius around the circle natal and settle locations
DispBuffer <- 20

### New Coordinate System
DispWide$natalX.new <- DispWide$natalX - min(trapcoords$UTMx) + BufferSize
DispWide$settleX.new <- DispWide$settleX - min(trapcoords$UTMx) + BufferSize
DispWide$natalY.new <- DispWide$natalY - min(trapcoords$UTMy) + BufferSize
DispWide$settleY.new <- DispWide$settleY - min(trapcoords$UTMy) + BufferSize
DispLong$Xnew <- DispLong$UTMx - min(trapcoords$UTMx) + BufferSize
DispLong$Ynew <- DispLong$UTMy - min(trapcoords$UTMy) + BufferSize

### Assign dispersers to Region 1 (south) or Region 2 (north) 
DispWide$natalRegion <- 1 + (DispWide$natalY >= 4262500)
DispWide$settleRegion <- 1 + (DispWide$settleY >= 4262500)
DispWide$Region <- DispWide$natalRegion * (DispWide$natalRegion == DispWide$settleRegion)


###############################################################################################################
### Plot traps & Dispersal Segments
###############################################################################################################
#ggplot(trapcoords, aes(x=Xnew, y=Ynew)) + 
#	geom_point(alpha=0.5, colour="red", shape=2, size=1) +
#	scale_shape(solid = FALSE) +
#	theme(panel.background=element_rect(fill="transparent", colour="NA")) +
#	coord_fixed() +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.min,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.min,y=S1y.max,xend=S1x.max,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.max,y=S1y.max,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.min,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.min,y=S2y.max,xend=S2x.max,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.max,y=S2y.max,xend=S2x.max,yend=S2y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.max,yend=S2y.min)) +
#	geom_line(data=DispLong, aes(x=Xnew,y=Ynew, group=ID), colour="blue")
#

###############################################################################
### Make sure Natal AND Settle points (with DispBuffer) are in S
### Removed (2013-12-19): AND make sure that the distance between the natal(X,Y) and settle(X,Y) is at least DispBuffer
###############################################################################

#DispWide <- subset(DispWide,	DispWide$distance > DispBuffer &
#				((	DispWide$natalX.new-DispBuffer > S1x.min & DispWide$natalX.new+DispBuffer < S1x.max &
#	   				DispWide$natalY.new-DispBuffer > S1y.min & DispWide$natalY.new+DispBuffer < S1y.max &
#					DispWide$settleX.new-DispBuffer > S1x.min & DispWide$settleX.new+DispBuffer < S1x.max &
#	   				DispWide$settleY.new-DispBuffer > S1y.min & DispWide$settleY.new+DispBuffer < S1y.max ) |
#				(	DispWide$natalX.new-DispBuffer > S2x.min & DispWide$natalX.new+DispBuffer < S2x.max &
#	   				DispWide$natalY.new-DispBuffer > S2y.min & DispWide$natalY.new+DispBuffer < S2y.max &
#					DispWide$settleX.new-DispBuffer > S2x.min & DispWide$settleX.new+DispBuffer < S2x.max &
#	   				DispWide$settleY.new-DispBuffer > S2y.min & DispWide$settleY.new+DispBuffer < S2y.max )))

DispWide <- subset(DispWide,	
				((	DispWide$natalX.new-DispBuffer > S1x.min & DispWide$natalX.new+DispBuffer < S1x.max &
	   				DispWide$natalY.new-DispBuffer > S1y.min & DispWide$natalY.new+DispBuffer < S1y.max &
					DispWide$settleX.new-DispBuffer > S1x.min & DispWide$settleX.new+DispBuffer < S1x.max &
	   				DispWide$settleY.new-DispBuffer > S1y.min & DispWide$settleY.new+DispBuffer < S1y.max ) |
				(	DispWide$natalX.new-DispBuffer > S2x.min & DispWide$natalX.new+DispBuffer < S2x.max &
	   				DispWide$natalY.new-DispBuffer > S2y.min & DispWide$natalY.new+DispBuffer < S2y.max &
					DispWide$settleX.new-DispBuffer > S2x.min & DispWide$settleX.new+DispBuffer < S2x.max &
	   				DispWide$settleY.new-DispBuffer > S2y.min & DispWide$settleY.new+DispBuffer < S2y.max )))

###############################################################################
### Subset Sessions (max was 36)
############################################################################### 

if(RunType == 1){
	# KeepSessions <- c(1:3,7:15)
	# KeepSessions <- c(7:15,20:24)
		# KeepSessions <- union(c(1:3,7:12,20:23), c(1,2,5:15,20,22:24))
		KeepSessions <- c(1:3, 5:15,20:24)	### revised 2014-07-06
}else if(RunType == 2){
	if(DisperserSet == 1){ 
		KeepSessions <- c(1:3,7:12,20:23)
	}else if(DisperserSet == 2){ 
		KeepSessions <- c(1,2,5:15,20,22:24)
	}else if(DisperserSet == 3){
		# KeepSessions <- union(c(1:3,7:12,20:23), c(1,2,5:15,20,22:24))
		KeepSessions <- c(1:3, 5:15,20:24)	### revised 2014-07-06
	}
}

### Subset
Hmatrix <- Hmatrix[Hmatrix$session %in% KeepSessions,]
SessionInfo <- SessionInfo[SessionInfo$sessionID %in% KeepSessions,]
DispWide <- subset(DispWide, natalsession %in% KeepSessions & settlesession %in% KeepSessions)

### Hold original session numbers
Hmatrix$session.original <- Hmatrix$session
SessionInfo$sessionID.original <- SessionInfo$sessionID
DispWide$natalsession.original <- DispWide$natalsession
DispWide$settlesession.original <- DispWide$settlesession

### Renumber individuals
Hmatrix$id.original <- Hmatrix$id
j <- 1
Hmatrix$id[1] <- 1
for(i in 2:nrow(Hmatrix)){
	if(Hmatrix$id.original[i] == Hmatrix$id.original[i-1]){
		Hmatrix$id[i] <- j
	}else{
		j <- j+1
		Hmatrix$id[i] <- j
	}
}



### Renumber Sessions
SessionInfo$sessionID <- 1:nrow(SessionInfo)
j <- 1
Hmatrix$session[1] <- 1
for(i in 2:nrow(Hmatrix)){
	if(Hmatrix$session.original[i] == Hmatrix$session.original[i-1]){
		Hmatrix$session[i] <- j
	}else{
		j <- j+1
		Hmatrix$session[i] <- j
	}
}

### Renumber sessions for Dispersal data
for(r in 1:nrow(DispWide)){
	for(s in 1:length(KeepSessions)){

		if(DispWide$natalsession.original[r] == KeepSessions[s]) { DispWide$natalsession[r] <- s }
		if(DispWide$settlesession.original[r] == KeepSessions[s]) { DispWide$settlesession[r] <- s }
	}
}


### Make New LONG Dispersal Data frame from the subsetted WIDE dataframe
DLr <- 2*nrow(DispWide)
DispLongSubset <- data.frame(Period=rep(NA,DLr), ID=rep(NA,DLr), Xnew=rep(NA,DLr), Ynew=rep(NA,DLr), session=rep(NA,DLr))
Lrow <- 0
for(r in 1:nrow(DispWide)){

	### Natal
		Lrow <- Lrow + 1

		DispLongSubset$Period[Lrow] <- "natal"
		DispLongSubset$ID[Lrow] <- DispWide$ID[r]
		DispLongSubset$Xnew[Lrow] <- DispWide$natalX.new[r]
		DispLongSubset$Ynew[Lrow] <- DispWide$natalY.new[r]
		DispLongSubset$session[Lrow] <- DispWide$natalsession[r]

	### Settle
		Lrow <- Lrow + 1

		DispLongSubset$Period[Lrow] <- "settle"
		DispLongSubset$ID[Lrow] <- DispWide$ID[r]
		DispLongSubset$Xnew[Lrow] <- DispWide$settleX.new[r]
		DispLongSubset$Ynew[Lrow] <- DispWide$settleY.new[r]
		DispLongSubset$session[Lrow] <- DispWide$settlesession[r]

}





###############################################################################################################
### Plot traps & Dispersal Segments -- Subset to KeepSessions
###############################################################################################################
#ggplot(trapcoords, aes(x=Xnew, y=Ynew)) + 
#	geom_point(alpha=0.5, colour="red", shape=2, size=1) +
#	scale_shape(solid = FALSE) +
#	theme(panel.background=element_rect(fill="transparent", colour="NA")) +
#	coord_fixed() +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.min,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.min,y=S1y.max,xend=S1x.max,yend=S1y.max)) +
#	geom_segment(aes(x=S1x.max,y=S1y.max,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S1x.min,y=S1y.min,xend=S1x.max,yend=S1y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.min,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.min,y=S2y.max,xend=S2x.max,yend=S2y.max)) +
#	geom_segment(aes(x=S2x.max,y=S2y.max,xend=S2x.max,yend=S2y.min)) +
#	geom_segment(aes(x=S2x.min,y=S2y.min,xend=S2x.max,yend=S2y.min)) +
#	geom_line(data=DispLongSubset, aes(x=Xnew,y=Ynew, group=ID), colour="blue")
#

###############################################################################
### Get an index of Session-Region
###############################################################################
Hmatrix$SessionRegion <- paste(Hmatrix$session,Hmatrix$group,sep="-")
Hmatrix$SR <- as.numeric(as.factor(Hmatrix$SessionRegion))


###############################################################################
### Get needed values from data and set other data values
###############################################################################

### Number of primary trapping sessions (36)
T.primary <- length(unique(Hmatrix$session))

### Number of Trapping Nights for each primary session
Trap.Nights <- SessionInfo$nights

### Number of traps
Num.Traps <- max(trapcoords$trapID)

### Number of Encounters
NumCap <- nrow(Hmatrix)

### Number of Individuals Caught 
### (same indiv caught in different primary session counts as new individual)
NumIndivsCaught <- max(Hmatrix$id)

### Number of individuals Caught by session
NumIndivsCaughtBySession <- rep(NA,T.primary)
for(s in 1:T.primary){
	NumIndivsCaughtBySession[s] <- length(unique(Hmatrix$id[Hmatrix$session==s]))
}

### Number of individuals Caught by Session-Region
NumIndivsCaughtBySessionRegion <- rep(NA,T.primary*2)
for(sr in 1:(2*T.primary)){
	NumIndivsCaughtBySessionRegion[sr] <- length(unique(Hmatrix$id[Hmatrix$SR==sr]))
}


### Number of individuals Caught by Region
NumIndivsCaughtByRegion <- rep(NA,2)
for(r in 1:2){
	NumIndivsCaughtByRegion[r] <- length(unique(Hmatrix$id[Hmatrix$group==r]))
}


###############################################################################################################
### Priors
### In email dated 2013-10-16, Karen Mabry: 
### "Just FYI, my very rough estimates for these sites during June-August 2005 are 57.5/ha for the "single" 
###	site, and 35 and 48.75/ha for the two grids included in the second analysis. "
###############################################################################################################
if(RunType == 1){
	M.density.ha <- 50	
}else if(RunType == 2){
	M.density.ha <- 100 
}
M1 <- round(M.density.ha * Area.S1.Hectares)
M2 <- round(M.density.ha * Area.S2.Hectares)


M <- rep(c(M1,M2), T.primary)

### need to have the index number of session and of Session-Region to go along with each "individual" 
###	(incl those not caught)

	### First get observed region, sessions for each indiv caught
	indiv.SesReg <- aggregate(cbind(session, group, SR) ~ id, data=Hmatrix, min)

	### Now tag on for those not caught
	M.SessionNum <- c(indiv.SesReg$session,   rep(rep(1:T.primary,each=2),	M-NumIndivsCaughtBySessionRegion)  )
	M.SesRegNum <-  c(indiv.SesReg$SR,		rep(1:(2*T.primary), 		M-NumIndivsCaughtBySessionRegion)  )
	M.RegionNum <-  c(indiv.SesReg$group,	rep(rep(1:2,T.primary),		M-NumIndivsCaughtBySessionRegion)  )




### Priors on Activity Centers
### For those caught in 1 trap, initialize to trap location + a random x and y
### For those caught in >1 trap, initialize to average trap location + a random x and y
### For those not caught, initialize to a uniform over S

### noise = random variation around initial activity centers for those indivs caught
noise  <- 4


SX.init <- array(NA, sum(M))
SY.init <- array(NA, sum(M))
n.caught <- array(0,sum(M))

### Fill for those that have been caught
H.row <- 0
for(i in 1:NumIndivsCaught){

	H.row <- H.row + 1
	x.temp <- trapcoords$Xnew[Hmatrix$trap[H.row]]
	y.temp <- trapcoords$Ynew[Hmatrix$trap[H.row]]
	n.caught.temp <- 1

	while(Hmatrix$id[H.row+1] == Hmatrix$id[H.row] & H.row < nrow(Hmatrix)){
		H.row <- H.row + 1
		n.caught.temp <- n.caught.temp + 1
		x.temp <- x.temp + trapcoords$Xnew[Hmatrix$trap[H.row]]
		y.temp <- y.temp + trapcoords$Ynew[Hmatrix$trap[H.row]]
	}

	n.caught[i] <- n.caught.temp
	SX.init[i] <- x.temp/n.caught.temp + runif(1, -1*noise, noise)
	SY.init[i] <- y.temp/n.caught.temp + runif(1, -1*noise, noise)

}



### Fill for those that have NOT been caught
SX.init[is.na(SX.init) & M.RegionNum==1] <- runif(n=(M1*T.primary-NumIndivsCaughtByRegion[1]), min=S1x.min, max=S1x.max)
SY.init[is.na(SY.init) & M.RegionNum==1] <- runif(n=(M1*T.primary-NumIndivsCaughtByRegion[1]), min=S1y.min, max=S1y.max)

SX.init[is.na(SX.init) & M.RegionNum==2] <- runif(n=(M2*T.primary-NumIndivsCaughtByRegion[2]), min=S2x.min, max=S2x.max)
SY.init[is.na(SY.init) & M.RegionNum==2] <- runif(n=(M2*T.primary-NumIndivsCaughtByRegion[2]), min=S2y.min, max=S2y.max)


###############################################################################################################
###############################################################################################################
### ANALYSIS
###############################################################################################################
###############################################################################################################

###############################################################################################################
### Model 
###############################################################################################################

### Specify Model in BUGS language
sink(paste(model.name,".txt",sep=""))
cat("
model{

	### Priors
		for(sr in 1:numSesReg){
			psi[sr] ~ dunif(0,1)
		}
		p0~dunif(0,1)
		sigma2~dunif(0,500)

	### Loop over all individuals (real and augmented)
	for(i in 1:M.total){

		### Loop over all traps
		for(j in 1:R){

			### D2 is the squared distance betw and activity center at (Sx,Sy) and trap j
			D2[i,j] <- pow(Sx[i]-TrapX[j],2) + pow(Sy[i]-TrapY[j],2) 

			K[i,j] <- exp(-1*D2[i,j] / sigma2)

			gamma[i,j] <- K[i,j]/E[i]
		}

		z[i] ~ dbern(psi[M.SesRegNum[i]])

		### Calculate the exposure for individual i
		E[i] <- sum(K[i,])

		### Number of captures for indiv i
		n[i] ~ dbin(Paug[i], T[M.SessionNum[i]])

		Paug[i] <- p[i]*z[i]

		p[i] <- p0*exp(1/-E[i])

		Sx[i] ~ dunif(xlower[M.RegNum[i]], xupper[M.RegNum[i]])
		Sy[i] ~ dunif(ylower[M.RegNum[i]], yupper[M.RegNum[i]])
	}

	for(h in 1:Ncap.total){
	
		### H1[h] is the trap number that indiv H2[h] encountered
		H1[h] ~ dcat(gamma[H2[h],])
	}



	### Derived Parameters - Natal & Settle Densities
	
		### d indexes the dispersers
		for(d in 1:numDispAll){
			for(i in 1:M.total){
				inNatal[d,i] <- z[i]*equals(natalsession[d],M.SessionNum[i])*((pow(Sx[i]-natalptsX[d],2)+pow(Sy[i]-natalptsY[d],2)) <= pow(DispBuffer,2))
				inSettle[d,i] <- z[i]*equals(settlesession[d],M.SessionNum[i])*((pow(Sx[i]-settleptsX[d],2)+pow(Sy[i]-settleptsY[d],2)) <= pow(DispBuffer,2))
	
				inNatalS[d,i] <- z[i]*equals(natalsession[d],M.SessionNum[i])
				inSettleS[d,i] <- z[i]*equals(settlesession[d],M.SessionNum[i])

			}
			### +1s needed because otherwise sometimes dividing by zero  +1 counts the disperser
			PctNChangeDisp[d] <- (sum(inSettle[d,1:M.total])+1) / (sum(inNatal[d,1:M.total]) +1)
			PctNChangeS[d] <- sum(inSettleS[d,1:M.total]) / sum(inNatalS[d,1:M.total]) 

			DensChangeDisp[d] <- (sum(inSettle[d,1:M.total]) - sum(inNatal[d,1:M.total])) / DispCircleAreaHectares
			DensChangeS[d] <-  (sum(inSettleS[d,1:M.total]) - sum(inNatalS[d,1:M.total]))  / RegionAreaHectares[DispRegion[d]]

			### 2014-06-13 add
			DensSettle[d] <- (sum(inSettle[d,1:M.total]) ) / DispCircleAreaHectares
			DensNatal[d]  <- (sum(inNatal[d,1:M.total]) ) / DispCircleAreaHectares
		}

		NatalN.1 <- sum(inNatal[1:numDisp1,1:M.total]) 
		SettleN.1 <- sum(inSettle[1:numDisp1,1:M.total]) 

		NatalNinS.1 <- sum(inNatalS[1:numDisp1,1:M.total]) 
		SettleNinS.1 <- sum(inSettleS[1:numDisp1,1:M.total]) 

		NatalN.2 <- sum(inNatal[(numDisp1+1):numDispAll,1:M.total]) 
		SettleN.2 <- sum(inSettle[(numDisp1+1):numDispAll,1:M.total]) 

		NatalNinS.2 <- sum(inNatalS[(numDisp1+1):numDispAll,1:M.total]) 
		SettleNinS.2 <- sum(inSettleS[(numDisp1+1):numDispAll,1:M.total]) 

		### Get OLS Regression Coeffs.  d indexes the dispersers
			### get parts 
			for(d in 1:numDispAll){
				xy.N[d] <- distance[d]*DensNatal[d]
			
				x2.N[d] <- DensNatal[d] * DensNatal[d]
			}

			### Slope for Disperser Set 1
			Beta.N1 <- (sum(xy.N[1:numDisp1]) - ((sum(DensNatal[1:numDisp1]) * sum(distance[1:numDisp1])) / numDisp1 )) / (  sum(x2.N[1:numDisp1]) - (( sum(DensNatal[1:numDisp1])*sum(DensNatal[1:numDisp1]) ) / numDisp1 ))

			### Slope for Disperser Set 2
			Beta.N2 <- (sum(xy.N[(numDisp1+1):numDispAll]) - ((sum(DensNatal[(numDisp1+1):numDispAll]) * sum(distance[(numDisp1+1):numDispAll])) / numDisp2 )) / (  sum(x2.N[(numDisp1+1):numDispAll]) - (( sum(DensNatal[(numDisp1+1):numDispAll])*sum(DensNatal[(numDisp1+1):numDispAll]) ) / numDisp2 ))

			### Slopes for each set WITHOUT the Dispersal Distance Outliers
			Beta.N1.NoOut <- (sum(xy.N[2:numDisp1]) - ((sum(DensNatal[2:numDisp1]) * sum(distance[2:numDisp1])) / (numDisp1-1) )) / (  sum(x2.N[2:numDisp1]) - (( sum(DensNatal[2:numDisp1])*sum(DensNatal[2:numDisp1]) ) / (numDisp1-1) ))



}
",fill = TRUE)
sink()



###############################################################################################################
### Prep/Package
###############################################################################################################


### Initial values
sigma2.post.max <- 300
inits <- function() list(psi=rep(1,(2*T.primary)), p0=runif(1,0,1), z=rep(1,sum(M)),
				Sx=SX.init, Sy=SY.init,
				sigma2=runif(1,0,sigma2.post.max) )

### Parameters Monitored
params <- c("p0", "sigma2", "psi", 
		"NatalN.1", "SettleN.1", "NatalNinS.1", "SettleNinS.1", 
		"NatalN.2", "SettleN.2", "NatalNinS.2", "SettleNinS.2", 
		"PctNChangeDisp", "PctNChangeS","DensChangeDisp", "DensChangeS", "DensSettle", "DensNatal", 
		"Beta.N1", "Beta.N2", "Beta.N1.NoOut" )

### data
DispCircleAreaHectares <- (pi * DispBuffer^2)/10000
numDisp1 <- length(DispWide$set[DispWide$set==1])
numDisp2 <- length(DispWide$set[DispWide$set==2])
data <- list(M.total=sum(M), Ncap.total=NumCap, TrapX=trapcoords$Xnew, TrapY=trapcoords$Ynew,
		H1=Hmatrix$trap, H2=Hmatrix$id, 
		M.SessionNum=M.SessionNum, M.RegNum=M.RegionNum, M.SesRegNum=M.SesRegNum,
		numSesReg=(2*T.primary), T=Trap.Nights, R=Num.Traps,
		xlower=c(S1x.min,S2x.min), xupper=c(S1x.max,S2x.max), 
		ylower=c(S1y.min,S2y.min), yupper=c(S1y.max,S2y.max),
		n=n.caught,
		natalptsX=DispWide$natalX.new, natalptsY=DispWide$natalY.new, natalsession=DispWide$natalsession,
		settleptsX=DispWide$settleX.new, settleptsY=DispWide$settleY.new, settlesession=DispWide$settlesession,
		numDisp1=numDisp1, numDisp2=numDisp2, numDispAll=nrow(DispWide),
		DispBuffer=DispBuffer, 
		DispCircleAreaHectares=DispCircleAreaHectares, 
		DispRegion=DispWide$Region, RegionAreaHectares=RegionAreaHectares,
		distance=DispWide$distance)

if(RunType == 1){
	### MCMC settings - trial
	ni=2000; nt=1; nb=300; nc=1
}else if(RunType == 2){
	### MCMC settings - actual
	ni=11000; nt=10; nb=1000; nc=3
}

###############################################################################################################
### Call JAGS from R
###############################################################################################################

	attempt <- 0

	print(Sys.time(), quote=FALSE)
	time.start <- Sys.time()

	out <- NULL

		while(is.null(out)){

			attempt <- attempt + 1
			cat("Attempt", attempt, "\tTime:",
				as.character(Sys.time()),"\n")
			flush.console()

			try(

				out <- jags(data, inits, params, paste(model.name,".txt",sep="")  ,
						n.chains=nc, n.thin=nt, n.iter=ni, n.burnin=nb)

			)
		}

	time.end <- Sys.time()
	print(Sys.time(), quote=FALSE)

	time.run <- time.end - time.start
	time.run

	### Save workspace
 	save.image(file.ws)




###############################################################################################################
### Results & Diagnostics
###############################################################################################################

print(out, dig=3)

### Average of (Change in Density around Disperser)-(Change in Density on Landscape)
AvgDiffChangeDensDispVsS.1 <- apply((out$BUGSoutput$sims.list$DensChangeDisp[,(1:numDisp1)]-out$BUGSoutput$sims.list$DensChangeS[,(1:numDisp1)]),1,mean)
AvgDiffChangeDensDispVsS.2 <- apply((out$BUGSoutput$sims.list$DensChangeDisp[,-(1:numDisp1)]-out$BUGSoutput$sims.list$DensChangeS[,-(1:numDisp1)]),1,mean)

### Posterior Distributions: p0, sigma2
# par(mfrow=c(2,1))
# hist(out$BUGSoutput$sims.list$p0, main="posterior p0",  xlab="", xlim=c(0,1))
# hist(out$BUGSoutput$sims.list$sigma2, main="posterior sigma2",  xlab="", xlim=c(0,sigma2.post.max))

### Posterior Distributions: Beta.N1, Beta.N2
# par(mfrow=c(2,1))
# hist(out$BUGSoutput$sims.list$Beta.N1 , main="posterior Beta.N1 ",  xlab="")
# hist(out$BUGSoutput$sims.list$Beta.N2, main="posterior Beta.N2",  xlab="")

### Posterior Distributions: Natal/Settle Density
# par(mfrow=c(2,2))
# rangeN1=range(c(out$BUGSoutput$sims.list$NatalN.1,out$BUGSoutput$sims.list$SettleN.1))
# rangeN2=range(c(out$BUGSoutput$sims.list$NatalN.2,out$BUGSoutput$sims.list$SettleN.2))
# hist(out$BUGSoutput$sims.list$NatalN.1 , main="posterior NatalN.1 ", xlim=rangeN1,  xlab="")
# hist(out$BUGSoutput$sims.list$NatalN.2 , main="posterior NatalN.2 ", xlim=rangeN2,  xlab="")
# hist(out$BUGSoutput$sims.list$SettleN.1, main="posterior SettleN.1",  xlim=rangeN1,  xlab="")
# hist(out$BUGSoutput$sims.list$SettleN.2, main="posterior SettleN.2",  xlim=rangeN2,  xlab="")

### Posterior Probability that 
mean(out$BUGSoutput$sims.list$NatalN.1 < out$BUGSoutput$sims.list$SettleN.1)
mean(out$BUGSoutput$sims.list$NatalN.1 > out$BUGSoutput$sims.list$SettleN.1)
mean(out$BUGSoutput$sims.list$NatalN.1 == out$BUGSoutput$sims.list$SettleN.1)
mean(out$BUGSoutput$sims.list$NatalN.2 < out$BUGSoutput$sims.list$SettleN.2)
mean(out$BUGSoutput$sims.list$NatalN.2 > out$BUGSoutput$sims.list$SettleN.2)
mean(out$BUGSoutput$sims.list$NatalN.2 == out$BUGSoutput$sims.list$SettleN.2)

### Change in N/density as ratio (Percent Change) 
### Posterior probs for percent change in N for dispersers vs. overall sample space
for(d in 1:nrow(DispWide)){
	print(mean(out$BUGSoutput$sims.list$PctNChangeDisp[,d] < out$BUGSoutput$sims.list$PctNChangeS[,d]))
}
### Posterior dist for relative change (ratio) in N for dispersers vs. overall sample space
for(d in 1:nrow(DispWide)){
	print(quantile((out$BUGSoutput$sims.list$PctNChangeDisp[,d] / out$BUGSoutput$sims.list$PctNChangeS[,d]),
		 probs=c(0.025, 0.5, 0.975)))
}


### Change in Density - Difference
### Posterior probs for percent change in N for dispersers vs. overall sample space
for(d in 1:nrow(DispWide)){
	print(mean(out$BUGSoutput$sims.list$DensChangeDisp[,d] < out$BUGSoutput$sims.list$DensChangeS[,d]))
}
### Posterior dist for relative change (ratio) in N for dispersers vs. overall sample space
for(d in 1:nrow(DispWide)){
	print(quantile((out$BUGSoutput$sims.list$DensChangeDisp[,d] - out$BUGSoutput$sims.list$DensChangeS[,d]),
		 probs=c(0.025, 0.5, 0.975)))
}


### Traceplot
# traceplot(out)


##################################################################################################
### Plot dispersal distance v. Natal Density
##################################################################################################

### First set-up data
DispWide$DensNatal.Median <-  apply(out$BUGSoutput$sims.list$DensNatal,2,median)
DispWide$DensNatal.025 <- apply(out$BUGSoutput$sims.list$DensNatal,2,quantile, probs = c(.025))
DispWide$DensNatal.975 <- apply(out$BUGSoutput$sims.list$DensNatal,2,quantile, probs = c(.975))

SetLabeller <- function(var, value){
	value <- as.character(value)
	if(var=="set"){
		value[value=="1"] <- "Telemetry"
		value[value=="2"] <- "Trapped"
	}
	return(value)
}

gg1 <- ggplot(DispWide, aes(x=DensNatal.Median,y=distance))+
	geom_point(alpha=0.3) +
	geom_errorbarh(aes(xmax = DensNatal.975, xmin = DensNatal.025), alpha = 0.3) +
	xlab("Local Natal Density (mice per ha)") +
	ylab("Dispersal Distance (m)") +
	theme(panel.background = element_blank(),
		panel.grid = element_blank(),
		panel.border=element_rect(fill=NA)) +
	facet_grid(set~.,  labeller=SetLabeller, scales="free_y") +
	stat_smooth(method="lm", se=FALSE, size=1.2) +
	expand_limits(y=0)

gg1.NoOut <- ggplot(DispWide[-1,], aes(x=DensNatal.Median,y=distance))+
	geom_point(alpha=0.3) +
	geom_errorbarh(aes(xmax = DensNatal.975, xmin = DensNatal.025), alpha = 0.3) +
	xlab("Local Natal Density (mice per ha)") +
	ylab("Dispersal Distance (m)") +
	theme(panel.background = element_blank(),
		panel.grid = element_blank(),
		panel.border=element_rect(fill=NA)) +
	facet_grid(set~.,  labeller=SetLabeller, scales="free_y") +
	stat_smooth(method="lm", se=FALSE, size=1.2) +
	expand_limits(y=0)
		





##################################################################################################
##################################################################################################
### Text output to text file
##################################################################################################
##################################################################################################
sink(outfile)
	cat(	"\nFilename:",outfile,
		"\n\n##################################################################################################
		### Input settings
		##################################################################################################",
		"\nDisperer Set =",DisperserSet,
		"\nBufferSize (meters) =",BufferSize,
		"\nArea.S.Hectares=",Area.S.Hectares,
		"\nArea.S1.Hectares=",Area.S1.Hectares,
		"\nArea.S2.Hectares=",Area.S2.Hectares,
		"\n\nDispBuffer (meters) =",DispBuffer,
		"\n\nKeepSessions =",KeepSessions,
		"\n\nPrior Density (per hectare), M.density.ha =",M.density.ha,
		"\nPrior on sigma =",sigma2.post.max,

		"\n\n##################################################################################################
		### Dispersers
		##################################################################################################\n"	
	)


	DispWide[,c("ID","distance","natalX",  "natalY", "natalsession",  "settlesession", "Region", "set")]

	cat(	"\n\n##################################################################################################
		### Summary output
		##################################################################################################\n"	
	)

	print(out, dig=3)

	cat(	"\n\n##################################################################################################
		### Posterior of Least Squares Slope y=Dispersal Distance, x=Natal Density 
		##################################################################################################\n"	
	)

	cat("\n### Set 1")
	cat("\n\tPosterior Prob Beta.N1 < 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N1 < 0)
	cat("\n\tPosterior Prob Beta.N1 > 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N1 > 0)

	cat("\n### Set 2")
	cat("\n\tPosterior Prob Beta.N2 < 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N2 < 0)
	cat("\n\tPosterior Prob Beta.N2 > 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N2 > 0)


	cat("\n\n### Set 1 - w/o Outliers")
	cat("\n\tPosterior Prob Beta.N1.NoOut < 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N1.NoOut < 0)
	cat("\n\tPosterior Prob Beta.N1.NoOut > 0\t\t")
	mean(out$BUGSoutput$sims.list$Beta.N1.NoOut > 0)



	cat(	"\n\n##################################################################################################
		### Posteriors Comparing Natal and Settle Densities - Over ALL Dispersers
		##################################################################################################\n"	
	)

	cat("\n### Set 1")
	cat("\n\tPosterior Prob NatalN.1 < SettleN.1\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.1 < out$BUGSoutput$sims.list$SettleN.1)
	cat("\tPosterior Prob NatalN.1 > SettleN.1\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.1 > out$BUGSoutput$sims.list$SettleN.1)
	cat("\tPosterior Prob NatalN.1 = SettleN.1\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.1 == out$BUGSoutput$sims.list$SettleN.1)

	cat("\n### Set 2")
	cat("\n\tPosterior Prob NatalN.2 < SettleN.2\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.2 < out$BUGSoutput$sims.list$SettleN.2)
	cat("\tPosterior Prob NatalN.2 > SettleN.2\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.2 > out$BUGSoutput$sims.list$SettleN.2)
	cat("\tPosterior Prob NatalN.2 = SettleN.2\t\t")
	mean(out$BUGSoutput$sims.list$NatalN.2 == out$BUGSoutput$sims.list$SettleN.2)

	cat(	"\n\n##################################################################################################
		### Change in N as ratio (Percent Change)
		##################################################################################################\n",
		"\nPosterior probs for percent change in N for dispersers vs. overall sample space\n"
	)
	for(d in 1:nrow(DispWide)){
		print(mean(out$BUGSoutput$sims.list$PctNChangeDisp[,d] < out$BUGSoutput$sims.list$PctNChangeS[,d]))
	}

	cat(	"\n\nPosterior Distribution for relative change (ratio) in N for dispersers vs. overall sample space\n"
	)
	for(d in 1:nrow(DispWide)){
	print(quantile((out$BUGSoutput$sims.list$PctNChangeDisp[,d] / out$BUGSoutput$sims.list$PctNChangeS[,d]),
		 probs=c(0.025, 0.5, 0.975)))
	}


	cat(	"\n\n##################################################################################################
		### Change in Density (# per Hectare) - Difference
		##################################################################################################\n",
		"\nPosterior probs for absolute change in N for dispersers vs. overall sample space\n"
	)

	for(d in 1:nrow(DispWide)){
		print(mean(out$BUGSoutput$sims.list$DensChangeDisp[,d] < out$BUGSoutput$sims.list$DensChangeS[,d]))
	}


	cat(	"\n\nPosterior dist for absolute change (difference) for dispersers vs. overall sample space\n"
	)
	for(d in 1:nrow(DispWide)){
		print(quantile((out$BUGSoutput$sims.list$DensChangeDisp[,d] - out$BUGSoutput$sims.list$DensChangeS[,d]),
			 probs=c(0.025, 0.5, 0.975)))
	}


	cat(	"\n\n##################################################################################################
		### Average of (Change in Density around Disperser)-(Change in Density on Landscape)
		### Disperser Set 1
		##################################################################################################\n"	)
		cat("\nSummary\n")
		summary(AvgDiffChangeDensDispVsS.1)

		cat("\n\n95% CrI:\t", quantile(AvgDiffChangeDensDispVsS.1,probs=0.025), "to", quantile(AvgDiffChangeDensDispVsS.1,probs=0.975))
		cat("\n\nProportion > 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.1 > 0)
		cat("\nProportion < 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.1 < 0)
		cat("\nProportion = 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.1 == 0)


	cat(	"\n\n##################################################################################################
		### Average of (Change in Density around Disperser)-(Change in Density on Landscape)
		### Disperser Set 2
		##################################################################################################\n"	)
		cat("\nSummary\n")
		summary(AvgDiffChangeDensDispVsS.2)

		cat("\n\n95% CrI:\t", quantile(AvgDiffChangeDensDispVsS.2,probs=0.025), "to", quantile(AvgDiffChangeDensDispVsS.2,probs=0.975))
		cat("\n\nProportion > 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.2 > 0)
		cat("\nProportion < 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.2 < 0)
		cat("\nProportion = 0:\t\t")
		mean(AvgDiffChangeDensDispVsS.2 == 0)



	cat(	"\n\n##################################################################################################
		### Run time
		##################################################################################################",
		"\nStart:\t", as.character(time.start),
		"\nEnd:\t", as.character(time.end),
		"\nRun Time (hours):\t", as.numeric(difftime(time.end, time.start, units="hours")),
		"\nRun Time (days):\t", as.numeric(difftime(time.end, time.start, units="days"))
	)
	
sink()



##################################################################################################
##################################################################################################
### Graphs output to PDF files
##################################################################################################
##################################################################################################
pdf(file=paste(model.name,"_",date.today,"_plots				.pdf", sep=""))


	### Posterior Distributions: p0, sigma2
	par(mfrow=c(2,1))
	hist(out$BUGSoutput$sims.list$p0, main="posterior p0",  xlab="", xlim=c(0,1))
	hist(out$BUGSoutput$sims.list$sigma2, main="posterior sigma2",  xlab="", xlim=c(0,sigma2.post.max))

	### Posterior Distributions: Beta.N1, Beta.N2
	par(mfrow=c(2,2))
	hist(out$BUGSoutput$sims.list$Beta.N1 , main="posterior Beta.N1 ",  xlab="")
	hist(out$BUGSoutput$sims.list$Beta.N2, main="posterior Beta.N2",  xlab="")
	hist(out$BUGSoutput$sims.list$Beta.N1.NoOut , main="posterior Beta.N1.NoOut ",  xlab="")


	### Posterior Distributions: Natal/Settle Density
	par(mfrow=c(2,2))
	rangeN1=range(c(out$BUGSoutput$sims.list$NatalN.1,out$BUGSoutput$sims.list$SettleN.1))
	rangeN2=range(c(out$BUGSoutput$sims.list$NatalN.2,out$BUGSoutput$sims.list$SettleN.2))
	hist(out$BUGSoutput$sims.list$NatalN.1 , main="posterior NatalN.1 ", xlim=rangeN1,  xlab="")
	hist(out$BUGSoutput$sims.list$NatalN.2 , main="posterior NatalN.2 ", xlim=rangeN2,  xlab="")
	hist(out$BUGSoutput$sims.list$SettleN.1, main="posterior SettleN.1",  xlim=rangeN1,  xlab="")
	hist(out$BUGSoutput$sims.list$SettleN.2, main="posterior SettleN.2",  xlim=rangeN2,  xlab="")


	gg1

	gg1.NoOut

dev.off()








##################################################################################################
##################################################################################################
### END OF PROGRAM - Save Workspace!
##################################################################################################
##################################################################################################
### Save workspace
 save.image(file.ws)