#####################################
#
# Author: Cody Evers
# Edited by: Gabriel Yospin
# Last revision date: 10.13.2011 0920 HRS PST
#
# This file will run simulations of probabilistic vegetation changes
# and then plot them. The number of simulations, years to simulate, and
# the parameters for the vegetation model are all user-specified.
#
#####################################

# Choose the appropriate directory
setwd(working_directory_name)

# Read in the necessary files
STM_master = read.csv(transition_table)
CTSS_aux = read.csv(transition_outcomes_table)
# Modify the ctss lookup to create a more condensed plot
CTSS_aux$yvalue = c(1:58)

# Specify the values of parameters for the model run
STM_sub = subset(STM_master, conserve == 0 
				& si == 100 
				& pvt == 6 
				& regen == 1)
STM_sub = STM_sub[,c("vegclass1","vegclass2","p")]
startState = 243 # Use 203 or 243
numRuns = 100
numYears = 500
list = c(1:(numRuns * numYears))
df = data.frame(matrix(list,nrow=numRuns,ncol=numYears))
colnames(df) = c(1:numYears)
df = df * NA
df[,1] = startState
for (i in c(1:numRuns)) {
	print(i)
	for (j in c(1:(numYears-1))) {
		veg1 = df[i,j]
		probTab =  subset(STM_sub, vegclass1 == veg1)
		mlist = NA
		for (x in c(1:nrow(probTab))) {
			list = rep(probTab$vegclass2[x], probTab$p[x] * 10000)
			mlist = c(mlist, list)
		}
		mlist = mlist[2:length(mlist)]
		veg2 = sample(mlist, 1)
		df[i,j+1] = veg2
	}
}
# Replace the veg class values with a condensed list from the lookup table
for (i in c(1:numYears)) {
	colData = data.frame(df[,i])
	colData = merge(colData, CTSS_aux, by.x = colnames(colData[1]), by.y = "vegstate1", all.x = TRUE)
	df[,i] = colData$yvalue
}
# Plot the results
plot(x = c(1,numYears), y = c(0,60), main = c("Successional Trajectories for ",startState," Over Time"),
     xlab = "Time (years)", ylab = "Vegetation Class", xlim = c(0,600), ylim = c(0,60), axes = FALSE, type="n")
# Color the background by veg type
rect(0, 0, 500, 6, border = NA, col=rgb(250,250,0,50,maxColorValue=255))#Yellow OA
rect(0, 6, 500, 10, border = NA, col=rgb(250,125,0,50,maxColorValue=255))#Orange OW
rect(0, 10, 500, 12, border = NA, col=rgb(180,100,0,50,maxColorValue=255))#Dk Orange OD
rect(0, 12, 500, 18, border = NA, col=rgb(140,250,0,50,maxColorValue=255))#Lt Green DO
rect(0, 18, 500, 27, border = NA, col=rgb(60,100,0,50,maxColorValue=255))#Dk Green DD
rect(0, 27, 500, 34, border = NA, col=rgb(50,175,130,50,maxColorValue=255))#Teal DM
rect(0, 34, 500, 36, border = NA, col=rgb(50,100,100,50,maxColorValue=255))#Slate DG
rect(0, 36, 500, 41, border = NA, col=rgb(0,250,250,50,maxColorValue=255))#Lt Blue BM
rect(0, 41, 500, 47, border = NA, col=rgb(250,180,250,50,maxColorValue=255))#Lt Pink PS
rect(0, 47, 500, 50, border = NA, col=rgb(80,10,60,50,maxColorValue=255))#Maroon PW
rect(0, 50, 500, 56, border = NA, col=rgb(200,80,250,50,maxColorValue=255))#Pink M
rect(0, 56, 500, 59, border = NA, col=rgb(200,0,100,50,maxColorValue=255))#Dk Pink MD
# Add an axis for time, with ticks every 100 years
axis(1, at = c(0,100,200,300,400,500), labels = c("0","100","200","300","400","500"),
     pos = 0.025, tck = 0.908, col.ticks = "white")
#axis(2, lwd = 0, lwd.ticks = 1, col.ticks = "grey", at = c(1:59), labels = NA, tck = 0.77, pos = 0.11)
## Add a legend
#leg.txt <- c("MD","M","PW","PS","BM","DG","DM","DD","DO","OD","OW","OA")
#legend("topright", legend = leg.txt, inset = 0, title = "Community Types",
#       fill = c(rgb(200,0,100,50,maxColorValue=255), rgb(200,80,250,50,maxColorValue=255), 
#		rgb(80,10,60,50,maxColorValue=255), rgb(250,180,250,50,maxColorValue=255),
#		rgb(0,250,250,50,maxColorValue=255), rgb(50,100,100,50,maxColorValue=255),
#		rgb(50,175,130,50,maxColorValue=255), rgb(60,100,0,50,maxColorValue=255),
#		rgb(140,250,0,50,maxColorValue=255), rgb(180,100,0,50,maxColorValue=255),
#		rgb(250,125,0,50,maxColorValue=255), rgb(250,250,0,50,maxColorValue=255)),
#      bty = "n")
# Add horizontal white lines
boxbounds = c(6,10,12,18,27,34,36,41,47,50,56)
for (i in c(1:length(boxbounds))) {
        abline(boxbounds[i], 0, col = "white", lwd = 3)				
}

# Plot the data
for (z in c(1:nrow(df))) {
	x = 1:(numYears)
	y = as.double(df[z,])
	lines(x = x, y = jitter(y, factor = 1, amount = 0.4),
	      col=rgb(0,0,0,10,maxColorValue=255),lwd=2)
}
# End