#######################################################################
## Nolte, E.G. and B. P. Pauli.  2011-2016  Raptor Migration Watch-site Simulations.  Code for R Stat. Envir. ##
#######################################################################

library(doParallel)
library(foreach)
setwd("/home")

#####Data import and format
Raptors <- read.delim("SS_Data.txt")		   ###covariate dataset for species of interest only, specify file here.
Raptors$Index <- c(1:length(Raptors$ID))
Model <- read.delim("TopModelBetas.txt")  ###specify table of coefficients
rownames(Model) <- Model$X
Model <- Model[,-1]
Plan <- read.delim("Plan.txt")				    ###a table of parameters to use in each simulation (each row specifies a number of parameters, and the number of iterations to run)

#####Set experiment via row in Plan
exp_num = 1

#####Specify desired directory for external files generated by this program, Do not include slash at end.
Output.dir <- '/home/output'

##functions to reconstitute p: change this code if model structure changes.  u= dataframe with covariates, v= beta value for observers.  
#reconstituteP is for 2 observers, used for Monte Carlo annual counts
#reconstituteP4 is for 4 observers (conditions in which the data was observed) used only to generate Horvitz-Thompson correction weights of available population.

model.constant <- Model['pInt','Beta']+Model['Accip','Beta']+(0.017241379*Model['size','Beta'])  
## the preceeding line is different for each species.  'Accip' and 0.017241379 for Sharp-shinned Hawk; 'Circus' and 0.362068966 for Northern Harrier.

reconstituteP <- function(r=Model, u, v, z=model.constant){
	indicator <- z+v+(sum(u[,c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2')]*r[c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2'),'Beta']))
	return((exp(indicator))/(1+exp(indicator)))
}
reconstituteP4 <- function(r=Model, u, z=model.constant, x){
	indicator1 <- z+max(Raptors[x, 2:11]*Model[2:11, 'Beta'])+(sum(u[,c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2')]*r[c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2'),'Beta']))
	indicator2 <- z+max(Raptors[x,12:21]*Model[2:11, 'Beta'])+(sum(u[,c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2')]*r[c('Alt', 'Alt_2', 'WindV', 'Day', 'Day_2'),'Beta']))
	p1 <- (exp(indicator1))/(1+exp(indicator1))
	p2 <- (exp(indicator2))/(1+exp(indicator2))
	return(p1+((1-p1)*p2))
}

##function to generate a year's count

yearrun <-  function(yf, nf=n, lf=l, kf=k, ef=e, ff=f, SampleSetf=Raptors$Index, bf=b, p.hat=Raptors$P4){	
	
	m <-ifelse(yf==1, nf, nf*(1+lf)^(yf-1))		#population change in year determined by specified trend
	mr <- m + rnorm(1, mean = 0, sd = m*kf)		#adds random residual variation
	Available <- SampleSetf[sample(x=1:length(SampleSetf), size= ifelse(test= mr>0, yes= mr, no=0), replace=TRUE, prob=1/p.hat)]
	
	P <- vector(length=length(Available), mode='numeric')
	
	if (length(Available)>=1) {
	
		#Calculate p values for year's observer effect
		for(x in c(1:length(Available))) {
		P[x] <- reconstituteP(u=Raptors[Available[x],], v=bf[yf])
		}
		
		#Monte Carlo proc to stochastically convert the individual probability of detection into a binomial outcome
		MonteCarlo <- vector(length=length(Available), mode= 'logical')
		for(h in c(1:length(P))){
			MonteCarlo[h] <- ifelse(runif(n=1, min=0, max=1) > P[h], yes=FALSE, no=TRUE)
		}
		return(c(length(Available),length(MonteCarlo[MonteCarlo==TRUE])))  #return the count
	} 
	else {return(c(0,0))
	}
}

##function to determine whether correct trend was detected 

trend.yes.no <- function(fyy=yy, fs=s, fYearCounts=YearCounts, C) {
	Trend <- lm(log(1+fYearCounts[1:fyy,C])~fYearCounts[1:fyy,"Year"])
	##extract confidence intervals as vectors (length=2)
	ci.Trend <- confint(Trend,level= fs)[2,]
	## Logical- determines whether confidence interval includes 0, collects in pre-defined vector. If true trend is a decline, use less than.  If increase, greater than.
	return(ifelse(test= (ci.Trend[1]<0 & ci.Trend[2]<0), yes=1, no=0))
}

##function to return trend coefficient
trend.coeff <- function(fyy=yy, fs=s, fYearCounts=YearCounts, C) {
	Trend <- lm(log(1+fYearCounts[1:fyy,C])~fYearCounts[1:fyy,"Year"])
	coeff <- coefficients(Trend)
	return(coeff[2])
}
	
P4 <- vector(length=nrow(Raptors), mode='numeric')   ##Pre-allocate memory for initial predictions of p
Raptors <- cbind(Raptors, P4)

registerDoParallel()   ##initate parallel backend- specify number of cores

##calculate estimated p for all observations
Raptors$P4 <- foreach(nn= 1:nrow(Raptors), .combine= "c", .inorder= TRUE) %dopar% {
	reconstituteP4(u=Raptors[nn,], x=nn)
}	

#####Experiment loop#####

for(j in exp_num:exp_num) {
	#####Set simulation parameters
	species <- "SSHA"	 				#string to identify species  
	i.set <- c(1:Plan[j,'ii'])   		#the number of iterations = 3000
	y.set <- c(1:Plan[j,'years'])	#the number of years per iteration = 30
	turnover <- Plan[j,"Turn"]		#probability of observer turnover. =1
	n <- Plan[j, "N"]					#starting population N1 = 100, 450, 2000
	l <- Plan[j, "pcy"]					#rate of population change per year = -0.035824
	s <- 1-Plan[j,"alpha"]			#1-alpha to detect trend = 0.9
	k <- Plan[j,"PopResid"]			#Coefficent of variability of population residuals CV(Nk)= 0.18, 0.26, 0.38
	e <- Plan[j,"Mean"]				#Observer norm dist mean
	f <- Plan[j,"SD"]					#Observer norm dist sd
	min.y <- 5							#minimum number of years to attempt trend analysis
	#####end section
	
	true <- foreach(y=y.set, .combine="c") %do% {
		ifelse(y==1, n, n*(1+l)^(y-1))
		}
	truetrend <- lm(log(1+true)~y.set)
	tt <- coefficients(truetrend)[2]
	rm(truetrend)
	
	tys <- foreach(o= min.y:max(y.set), .combine="c") %do% {				##makes column names for output frame
		paste("year",o, sep=".")
		}
		
	C.Trend.correct <- matrix(nrow=length(i.set), ncol=length(y.set)-(min.y-1), dimnames=list(NULL, tys)) 
	A.Trend.correct <- matrix(nrow=length(i.set), ncol=length(y.set)-(min.y-1), dimnames=list(NULL, tys)) 
	C.Trend.correct <- data.frame(C.Trend.correct)
	A.Trend.correct <- data.frame(A.Trend.correct)						###define output frame for analysis at all years 9/22/11 EN
	A.Coeff <- vector(length=length(i.set), mode="numeric")
	C.Coeff <- vector(length=length(i.set), mode="numeric")
	A.Coeff.all <- matrix(nrow=length(i.set), ncol=length(y.set)-(min.y-1), dimnames=list(NULL, tys))
	C.Coeff.all <- matrix(nrow=length(i.set), ncol=length(y.set)-(min.y-1), dimnames=list(NULL, tys))
	A.Coeff.all <- data.frame(A.Coeff.all)
	C.Coeff.all <- data.frame(C.Coeff.all)
	b <- vector(length=length(y.set), mode= 'numeric')					###define observer effect vector
	
	####Iteration Loop####
	for(i in i.set) {
		##select year observer effect value
		for(y in y.set){
			d <- ifelse(test= runif(n=1, min=0, max=1)<turnover, TRUE, FALSE)
			b[y] <- ifelse(test= y==1 | d==TRUE, yes= rnorm(n=1, mean= e, sd= f), no= b[(y-1)])
		}
	
		####Year Loop
		YearCounts <-  foreach(y= y.set, .combine= "rbind", .inorder= TRUE) %dopar% { 
			yearrun(yf=y)
		}
				
		Year <- vector(length=length(y.set), mode='integer')+y.set
		Iter <- vector(length=length(y.set), mode='integer')+i
		YearCounts <- data.frame(Avail=YearCounts[,1], Count=YearCounts[,2], Year, Iter, row.names=NULL)
		YearCounts[is.na.data.frame(YearCounts)] <- 0
	
		##Run linear regression
		C.Trend.correct [i,] <- foreach(yy= min.y:max(y.set), .inorder= TRUE) %dopar% {
			trend.yes.no(C="Count", fyy=yy)	
		}		
		A.Trend.correct[i,] <- foreach(yy= min.y:max(y.set), .inorder= TRUE) %dopar% {
			trend.yes.no(C="Avail", fyy=yy)	
		}
		A.Coeff.all[i,] <- foreach(yy= min.y:max(y.set), .inorder= TRUE) %dopar% {
			trend.coeff(C="Avail", fyy=yy) 	
		}
		C.Coeff.all[i,] <- foreach(yy= min.y:max(y.set), .inorder= TRUE) %dopar% {
			trend.coeff(C="Count", fyy=yy) 	
		}		
	print(paste("iteration",i,"complete"))
	write.table(C.Coeff.all[i,],file= paste(Output.dir, "/count_output_exp_", j, ".txt", sep=""), append=TRUE,  row.names=FALSE, col.names= ifelse(i==1, TRUE, FALSE),sep="\t")
	write.table(A.Coeff.all[i,],file= paste(Output.dir, "/avail_output_exp_", j, ".txt", sep=""), append=TRUE,  row.names=FALSE, col.names= ifelse(i==1, TRUE, FALSE),sep="\t")
	write.table(tt,file= paste(Output.dir, "/true_output_exp_", j, ".txt", sep=""), append=TRUE,  row.names=FALSE, col.names= ifelse(i==1, TRUE, FALSE),sep="\t")
	}
	
	###Output.
	Sp <- vector(length= max(y.set)-(min.y-1), mode='character') ; Sp[1:length(Sp)] <- species
	N <- vector(length= max(y.set)-(min.y-1), mode='integer')+n
	iterations <- vector(length= max(y.set)-(min.y-1), mode='integer')+length(i.set)
	min.years <- vector(length= max(y.set)-(min.y-1), mode='numeric')+min.y
	max.years <- vector(length= max(y.set)-(min.y-1), mode='integer')+length(y.set)
	PopResid <- vector(length= max(y.set)-(min.y-1), mode='numeric')+k
	pcy <- vector(length= max(y.set)-(min.y-1), mode='numeric')+l
	Turn <- vector(length= max(y.set)-(min.y-1), mode='numeric')+turnover
	Obs_Mean <- vector(length= max(y.set)-(min.y-1), mode='numeric')+e
	Obs_SD <- vector(length= max(y.set)-(min.y-1), mode='numeric')+f
	alpha <- vector(length= max(y.set)-(min.y-1), mode='numeric')+(1-s)
	A.Bias <- vector(length= max(y.set)-(min.y-1), mode='numeric')+mean(A.Coeff-tt)
	P.Bias <- vector(length= max(y.set)-(min.y-1), mode='numeric')+mean(C.Coeff-tt)
	Years <- vector(length= max(y.set)-(min.y-1), mode='integer')+c(min.y:max(y.set))
	Avail.Power <- foreach(z= 1:(max(y.set)-(min.y-1)), .combine="c") %dopar% {
			sum(A.Trend.correct[,z])/length(i.set)
		}
	Count.Power <- 	foreach(z= 1:(max(y.set)-(min.y-1)), .combine="c") %dopar% {
			sum(C.Trend.correct[,z])/length(i.set)
		}
	Output <- data.frame(Sp, N, iterations, min.years, max.years, PopResid, pcy, Turn, Obs_Mean, Obs_SD, alpha, A.Bias, P.Bias, Years, Avail.Power, Count.Power)
	write.table(Output, file= paste(Output.dir, "/output_exp_", j, ".txt", sep=""), append=TRUE,  row.names=FALSE, sep="\t")
}