# Scott Merrill
# Oct 12
# distribution of medical supplies

rm(list = ls(all=TRUE))


# Schistosomiasis is a parasitic disease propogated through contaminated water 
# Snails are an intermediary host of the parasite
# 30 advanced aid stations were arranged evenly spaced along a river
# Site 1 is at the river source, site 30 is at the river mouth
# The parasite is known to travel in open water
# Once detected, there is an approximately 70% chance that another outbreak will occur downriver
#	with a chance of occurrence following a beta distribution (skewed) 
#	with shape parameters 6 and 2.6 
# Secondary outbreaks occur at a distance away from the primary outbreak station 
#	with distance between primary and secondary outbreaks following a poisson distribution
#	with a lambda = 3 (mean and standard deviation = 3)

# about 3 outbreaks occur each year(poisson with a lambda of 3) along the river
# snail populations are more highly concentrated towards the the river source 
# the location of the initial outbreaks appear to also follow a poisson distribution 
#	with a lambda of 5 (site number)


# if resources are limited. Which 7 aid stations should receive treatment supplies? 

# goal: simulate system 10000 times to figure out which stations are most likely to observe outbreaks
# I will assume that an aid station can either observe an outbreak or not (0 or 1)
#	Thus, in any one simulation one station cannot observe more than one outbreak
#	If an outbreak occurs at a site that already has disease, it does not add to the resources needed


# initial disease matrix year 1

# to develop a probability of occurance at any station
# 30 river sites, 2 outbreak types, 10000 simulations
disease = array(data = 0, dim = c(30,2,10000))


# number of simulations
for (sim in 1:10000) {
	# Need to create initial outbreaks	
	number.outbreaks = rpois(lambda = 3,n=1)
	print(c("number of outbreaks = ",number.outbreaks))
	# place the outbreaks at sites
	primary.outbreaks = rpois(lambda = 5, n = number.outbreaks)
	print(c("location of outbreaks = ",primary.outbreaks))
	
	# place outbreaks into disease matrix	
	disease[primary.outbreaks,1,sim] = 1

	# create secondary outbreaks
	for (site in 1:30) {
		if (disease[site,1,sim] ==1) {
			occur=rbeta(n=1,shape1 = 6, shape2 = 2.6)
			if (occur < .7) {
				s.outbreak.location = rpois(lambda=3,n=1) + site
				# place outbreak unless it is past the river mouth (30th station)
				if (site < 31) {
					disease[s.outbreak.location,2,sim] = 1
				}
			}
		}
	}
	# making sure that I am not double counting 
	#	If a secondary outbreak occurs at where an outbreak is already occurring
	#	there should be only one outbreak. I am zeroing out my secondary outbreaks in this simulation.
	for (site in 1:30) {
		if (disease[site,1,sim] == 1 &  disease[site,2,sim] == 1) {
			disease[site,2,sim] = 0
		}
	}
}



disease[,,c(1:10)]
sum.incidence1=numeric(30)
total.incidence=numeric(30)
sum.incidence2=numeric(30)


for (site in 1:30) { 
	sum.incidence1[site]=sum(disease[site,1,c(1:10000)])
	sum.incidence2[site]=sum(disease[site,2,c(1:10000)])
	total.incidence[site]=sum(disease[site,1,c(1:10000)])+ sum(disease[site,2,c(1:10000)])
}

total.incidence

plot(sum.incidence1+sum.incidence2, type = "b", col = "red",pch = 1,lwd = 3, 
	main = "Bilharzia Outbreaks by Station")
points(sum.incidence1,col = "green", lwd = 3, cex = 1, type = "b")
points(sum.incidence2,col = "blue", lwd = 3, cex = 1, type = "b")

rect(22,3000,26.9,2445, col="bisque")
leg.txt <- c("Total", "Primary", "Secondary")
legend(x = 22, y = 3000,legend = leg.txt, col=c("red","green","blue"), cex =1, pch=1, lwd = 2,
         lty=1, merge=TRUE)#, trace=TRUE)

# Based on this simulation, I might recommend putting my 7 stashes of treatment supplies at stations 2-8