# Code by Matt Asher for statisticsblog.com
# Feel free to modify and redistribute, but please keep this notice 

costPerKafon = 50
	
# How many plates per kafon (it dies when hits this many)
plates = 4

# Changing this changes the variability in the Kafon's path
windSpeed = 1

# Dimmensions of mine field
depth = 1000
width = 1000

# Chance that a given location will have a mine
# We will try out diffent mine probabilities
# This may take a long time to run!
# To reduce run time, lower the number of mineProbs and trailsPerP
mineProbs = (1/2)^(10:23)

# How many trials should we simulate per mineProbs?
trialsPerP = 200

# Max cost we are willing to spend to check if a field has mines
costThreshold = 30000
maxKafons = floor(costThreshold/costPerKafon)

# Where will we be starting the kafons?
# Function to find optimal starting points when total number of kafons is 
# indeterminate. This may not be an efficient algorithm
buildStartPoints <- function(n, width) {
	startPoints = rep(0, n)
	k = 1
	i = 1
	count = 0
	while(k < n) {
		countBy = (1/2)^i
		count = countBy
		while(count < 1) {
			if(!count %in% startPoints) {
				startPoints[k] = count
				k = k + 1
			}
			count = count + countBy
		}
		i = i+1
			
	}
	# return(startPoints)
	
	return(round(startPoints*width))
}

startPoints = buildStartPoints(maxKafons, width)

# Create a blank data frame with trialsPerP * length(mineProbs) rows and 
# Columns for kafonsRun, minesFound, minesThere, 
results = matrix(0, nrow=(trialsPerP * length(mineProbs)), ncol=5)
colnames(results) = c("mineP", "kafonsRun", "mineFound", "minesThere", "areaCovered")

curRow = 0
for(mineP in mineProbs) {
	
	for(i in 1:trialsPerP) {
		# Generate the minefield
		covered = matrix(0, nrow=width, ncol=depth)
		mineMatrix = matrix(runif(width*depth), nrow=width, ncol=depth)
		mineMatrix = ifelse(mineMatrix < mineP, 1, 0)
		
		curRow = curRow + 1
		kafonsTried = 1
		mineFound = FALSE
		
		# Go until we find a mine or hit max cost
		while(!mineFound & (kafonsTried < maxKafons)) {
			startY = startPoints[i]
			startX = 0
			
			xPos = startX
			yPos = startY
		
			for(j in 1:depth) {
				
				# This is the noise in movement, feel free to adjust the function
				xPos = xPos + windSpeed
				yPos = yPos + rnorm(1,0,1)
				rndY = round(yPos)
				
				# Are we in the field?
				if(xPos > 0 && xPos <= depth && rndY > 0 && rndY <= width) {
					# Set this part of covered to 1
					if(covered[rndY,xPos] != 1) {
						covered[rndY,xPos] = 1
					}
				
					# Did we hit a mine? 
					if(mineMatrix[rndY,xPos] == 1) {
						
						# Note that we found a mine this trial
						results[curRow,] = c(mineP, kafonsTried, 1, sum(mineMatrix), sum(covered)) 
						
						# Note that at least one mine has been found
						mineFound = TRUE
						
						break;
					}
				}
			}
			kafonsTried = kafonsTried + 1
		}
		
		# If none found, record the data here
		if(!mineFound) {
			results[curRow,] = c(mineP, kafonsTried, 0, sum(mineMatrix), sum(covered))
		} 	
	}
}
 par(mar=c(5.1,4.1,4.1,1.0)) plot(.5,.5,ylim=c(0, costThreshold), 
 xlim=rev(c(min(mineProbs),max(mineProbs))), yaxt="n", col="white", log="x", 
 xlab="Mine probability per unit", ylab="Cost to detect first mine", bty="n") 
 title(main="Cost to detect the presence of landmines \n vs. probability of a 
 mine in a given zone. \n Numbers represent the true number of mines in the 
 field.", cex.main=1) axis(2,at=axTicks(2), labels=sprintf("$%s", axTicks(2))) # 
 points(results[1:2800,"mineP"], (results[1:2800,"kafonsRun"]*costPerKafon), 
 pch=20, col="blue")