# Scott Merrill
# Code for October 3rd lecture
# Sept 28, 2012
# Quantitative thinking in the life sciences

rm(list = ls(all = TRUE))


# Setting up the initial conditions for the array
idaho =array(0, dim= c(6,5,10))

# Initial tree density for each class
idaho[,4,1] = c(3,1,4,0,2,1) # mature trees in plots in year 1
idaho[,2,1]=c(1,0,1,0,0,1) # 1 year old saplings in plots in year 1
idaho[,3,1]=c(1,0,0,1,0,1) # 2 year old saplings in plots in year 1
idaho[,1,1]=c(0,0,0,0,0,0) # sprouts in plots in year 1

idaho[,,c(1:2)]


# sites that burned in year 1 receive a 1 in the fire column (column 5).
idaho[,5,1]=c(0,0,0,1,1,1) # fire in plots 4 to 6 in year 1
idaho[,5,4]=c(1,1,1,0,0,0) # fire in plots 1 to 3 in year 4
idaho[,5,7]=c(1,1,1,1,1,1) # fire in all plots in year 7



# adding sprouts because having no sprouts was no fun
####
#Note: These sprouts were not included in the assignment set-up
idaho[,1,1]=c(3,2,1,0,3,4) # sprouts in plots in year 1

idaho[,,c(1:2)]

# curiousity variable inserted to see how many trees are being killed
trees.killed = 0

# setting up my year loop. My tenth year is informed by my 9th year
#	so, I only have my time loop running from year 1 to year 9
for (t in 1:9) {
	# Setting up my site/plot loop. Running through all 6 plots
	for (y in 1:6) {
		#####
		# Setting up my fire if then statement loop
		#####
		if (idaho[y,5,t] == 1) {
			# kill sprouts and saplings if there is a fire
			idaho[y,c(1:3),t] = 0
			# sprouts emerge the following year, one for every mature tree
			idaho[y,1,t+1] = idaho[y,4,t]
		}
		#####
		# Starting growth code
		#####
		idaho[y,2,t+1] = idaho[y,1,t]
		idaho[y,3,t+1] = idaho[y,2,t]
		idaho[y,4,t+1] = idaho[y,3,t] + idaho[y,4,t]
	} # end y loop (site loop)
} # end t loop (time loop)

# look to see if example is working
idaho

###### 
#filling a matrix with values of the number of organisms in each plot across the years
# 	then plotting that matrix
######

org.time.matrix=matrix(nrow=10,ncol=6,data=0)

for (y in 1:6) {
	for (t in 1:10) { 
		org.time.matrix[t,y] = sum(idaho[y,c(1:4),t])
	} # end time fill loop
} # end site fill loop

matplot(org.time.matrix,type = "b")		



####################################################################################################
# restarting the exercise with beetle kill
##############################

rm(list = ls(all = TRUE))


# Setting up the initial conditions for the array
idaho =array(0, dim= c(6,5,10))

# Initial tree density for each class
idaho[,4,1] = c(3,1,4,0,2,1) # mature trees in plots in year 1
idaho[,2,1]=c(1,0,1,0,0,1) # 1 year old saplings in plots in year 1
idaho[,3,1]=c(1,0,0,1,0,1) # 2 year old saplings in plots in year 1
idaho[,1,1]=c(0,0,0,0,0,0) # sprouts in plots in year 1

# check my work
idaho[,,c(1:2)]

# sites that burned in year 1 receive a 1 in the fire column (column 5).
idaho[,5,1]=c(0,0,0,1,1,1) # fire in plots 4 to 6 in year 1
idaho[,5,4]=c(1,1,1,0,0,0) # fire in plots 1 to 3 in year 4
idaho[,5,7]=c(1,1,1,1,1,1) # fire in all plots in year 7

# adding sprouts because having no sprouts was no fun
####
#Note: These sprouts were not included in the assignment set-up
idaho[,1,1]=c(3,2,1,0,3,4) # sprouts in plots in year 1

# check my work
idaho[,,c(1:2)]

# curiousity variable inserted to see how many trees are being killed
trees.killed = 0

# setting up my year loop. My tenth year is informed by my 9th year
#	so, I only have my time loop running from year 1 to year 9
for (t in 1:9) {
	# Setting up my site/plot loop. Running through all 6 plots
	for (y in 1:6) {
		#####
		# Setting up my fire if then statement loop
		#####
		if (idaho[y,5,t] == 1) {
			# kill sprouts and saplings if there is a fire
			idaho[y,c(1:3),t] = 0
			# sprouts emerge the following year, one for every mature tree
			idaho[y,1,t+1] = idaho[y,4,t]
		}

		#####
		# Beetle kill loops
		#####
		# I am running my beetle kill loop after fire because I figure that mature trees
		#	have put out cones in the same year that they were killed by beetles, 
		#	allowing for sprouts
		# if statement to start checking if there are mature trees
		#	warning: a for loop reading: for (t in 1:0) runs explanatory code at bottom		
		if (idaho[y,4,t] > 0) {
			# create a variable that equals the number of mature trees in idaho[y,4,t]
			# this variable will be used to give my loop a number of times to cycle
			#	That is, each mature tree has a chance to be killed, in my loop, by the beetles
			# 	before beetles have had a chance to kill (and thus reduce) the number of trees
			mature = idaho[y,4,t]
			for (beetle in 1:mature) {
				# selecting a random integer,named kill, from 1 to rate
				rate = 10 # kill rate = 1/rate, in this case one in 10 trees is likely to die			
				kill = sample(1:rate,1)
				# creating an if statement that says that if kill = 3 
				#	then I reduce the number of mature trees by one this year (t)
				if (kill == 3) {
					# I added some print statements in my loop to see how many trees were dying
					#	and in which plots, at what time. 
					print(c("beetle kill in plot ",y," and year ",t))
					print(c("starting mature trees",idaho[y,4,t]))
					# given that a mature tree is killed, line below reduces the number of mature trees by 1
					idaho[y,4,t] = idaho[y,4,t] - 1
					print(c("ending mature trees",idaho[y,4,t]))
					trees.killed = trees.killed + 1
					print(c("trees.killed =",trees.killed))
					print("")
					
				} # ending "kill" if loop
			} # ending beetle for loop
		} # Ending idaho[y,4,t] > 0 if loop (i.e., do mature trees exist)
		#####
		# Starting growth code
		#####
		idaho[y,2,t+1] = idaho[y,1,t]
		idaho[y,3,t+1] = idaho[y,2,t]
		idaho[y,4,t+1] = idaho[y,3,t] + idaho[y,4,t]
	} # end y loop (site loop)
} # end t loop (time loop)

# look to see if example is working
idaho

###### 
#filling a matrix with values of the number of organisms in each plot across the years
# 	then plotting that matrix
######

org.time.matrix=matrix(nrow=10,ncol=6,data=0)

for (y in 1:6) {
	for (t in 1:10) { 
		org.time.matrix[t,y] = sum(idaho[y,c(1:4),t])
	} # end time fill loop
} # end site fill loop

matplot(org.time.matrix,type = "b")		


##################################################################################################
# same exercise but this time with random fires 
# fires occur at a rate of  1 in 5 chance of burning per year
# that is, each plot has a 20% chance of burning each year

rm(list = ls(all = TRUE))


# Setting up the initial conditions for the array
idaho =array(0, dim= c(6,5,10))

# Initial tree density for each class
idaho[,4,1] = c(3,1,4,0,2,1) # mature trees in plots in year 1
idaho[,2,1]=c(1,0,1,0,0,1) # 1 year old saplings in plots in year 1
idaho[,3,1]=c(1,0,0,1,0,1) # 2 year old saplings in plots in year 1
idaho[,1,1]=c(0,0,0,0,0,0) # sprouts in plots in year 1

idaho[,,c(1:2)]

##### 
# create fires for the entire array of sites and years
#####
# e.g., sites that randomly burned in year 1 receive a 1 in the fire column (column 5).

# establishing the fire rate
fire.rate = 5

# starting my random fire placement loop
for (y in 1:6) {
	for (t in 1:10) { 
		fire = sample(1:fire.rate,1)
		if (fire == 5) {
			idaho[y,5,t] = 1
		}
	}
}

idaho[,5,] # checking to see if it worked
sum(idaho[,5,]) # the sum should equal approximately 60/fire.rate (with a binomial distribution)



# adding sprouts because having no sprouts was no fun
####
#Note: These sprouts were not included in the assignment set-up
idaho[,1,1]=c(3,2,1,0,3,4) # sprouts in plots in year 1

idaho[,,c(1:2)]

# curiousity variable inserted to see how many trees are being killed
trees.killed = 0

# setting up my year loop. My tenth year is informed by my 9th year
#	so, I only have my time loop running from year 1 to year 9
for (t in 1:9) {
	# Setting up my site/plot loop. Running through all 6 plots
	for (y in 1:6) {
		#####
		# Setting up my fire if then statement loop
		#####
		if (idaho[y,5,t] == 1) {
			# kill sprouts and saplings if there is a fire
			idaho[y,c(1:3),t] = 0
			# sprouts emerge the following year, one for every mature tree
			idaho[y,1,t+1] = idaho[y,4,t]
		}

		#####
		# Beetle kill loops
		#####
		# I am running my beetle kill loop after fire because I figure that mature trees
		#	have put out cones in the same year that they were killed by beetles, 
		#	allowing for sprouts
		# if statement to start checking if there are mature trees
		#	warning: a for loop reading: for (t in 1:0) runs explanatory code at bottom		
		if (idaho[y,4,t] > 0) {
			# create a variable that equals the number of mature trees in idaho[y,4,t]
			# this variable will be used to give my loop a number of times to cycle
			#	That is, each mature tree has a chance to be killed, in my loop, by the beetles
			# 	before beetles have had a chance to kill (and thus reduce) the number of trees
			mature = idaho[y,4,t]
			for (beetle in 1:mature) {
				# selecting a random integer,named kill, from 1 to rate
				rate = 10 # kill rate = 1/rate, in this case one in 10 trees is likely to die			
				kill = sample(1:rate,1)
				# creating an if statement that says that if kill = 3 
				#	then I reduce the number of mature trees by one this year (t)
				if (kill == 3) {
					# I added some print statements in my loop to see how many trees were dying
					#	and in which plots, at what time. 
					print(c("beetle kill in plot ",y," and year ",t))
					print(c("starting mature trees",idaho[y,4,t]))
					# given that a mature tree is killed, line below reduces the number of mature trees by 1
					idaho[y,4,t] = idaho[y,4,t] - 1
					print(c("ending mature trees",idaho[y,4,t]))
					trees.killed = trees.killed + 1
					print(c("trees.killed =",trees.killed))
					print("")
					
				} # ending "kill" if loop
			} # ending beetle for loop
		} # Ending idaho[y,4,t] > 0 if loop (i.e., do mature trees exist)
		#####
		# Starting growth code
		#####
		idaho[y,2,t+1] = idaho[y,1,t]
		idaho[y,3,t+1] = idaho[y,2,t]
		idaho[y,4,t+1] = idaho[y,3,t] + idaho[y,4,t]
	} # end y loop (site loop)
} # end t loop (time loop)

# look to see if example is working
idaho

###### 
#filling a matrix with values of the number of organisms in each plot across the years
# 	then plotting that matrix
######

org.time.matrix=matrix(nrow=10,ncol=6,data=0)

for (y in 1:6) {
	for (t in 1:10) { 
		org.time.matrix[t,y] = sum(idaho[y,c(1:4),t])
	} # end time fill loop
} # end site fill loop

matplot(org.time.matrix,type = "b")		

##########################################################################################################
#### example code for Tom's population growth rate matrix and matplot
lam=seq(1,1.6,by=0.1)

data=matrix(nrow=10,ncol=length(lam),data=10)

for (time in 2:10) { 
	for (lambda in 1:length(lam)) { 
		data[time,lambda] = data[time-1,lambda]*lam[lambda]
	}
}
matplot(data,type = "b")




#########################################################################################
# r warning regarding for loops
###########
# compile the following


for (x in 1:0) {	
	print (x)
}
# so if you run a beetle kill loop, 
#	initialized by the number of mature trees in the cell, 
#	given no mature trees exist (i.e., the cell = 0),
# 	R will run through the loop twice! 
#	in my code that would potentially kill off two trees
#	if I had not used the if statement before starting my for loop

# compile this loop as well. R is just running through the most likely sequence given your input
for (x in 1:-2) {	
	print (x)
}