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

# set.seed(345)

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 = 500

# Chance that a given location will have a mine
mineP = 0.003

# We will try out Kafons from min to max
minKafons = 1
maxKafons = 2000

# Jump up, so you don't do every one in that range
toCheck = seq(minKafons, maxKafons, 20)

# Track variables over the trials
history.numbKafons = rep(0, length(toCheck))
history.minesLeft = rep(0, length(toCheck))
history.totalAreaCovered = rep(0, length(toCheck))
history.numbMines = rep(0,  length(toCheck))
# Number of kafons that use up all their plates before reaching the end of the minefield
history.platedOut = rep(0,  length(toCheck))

# Should we show the plotting? Only do this if you are trying out a single number of Kafons!
plotKafons = F

trial = 1
for(k in toCheck) {

	kafons = k
	
	if(plotKafons) {
		plot(0, 0, col="blue", pch=".", xlim=c(0,depth), ylim=c(0,width), xlab="Minefield depth", ylab="Minefield width", axes=F)
		Axis(side=1, at=c(1,depth), labels=c(1,depth))
		Axis(side=2, at=c(1,width), labels=c(1,width))
	}
	
	covered = matrix(0, nrow=width, ncol=depth)
	
	mineMatrix = matrix(runif(width*depth), nrow=width, ncol=depth)
	mineMatrix = ifelse(mineMatrix < mineP, 1, 0)
	
	startMines = sum(mineMatrix)
	
	# Track the number of kafons which stopped working while still in the field
	platedOut = 0
	
	for(i in 1:kafons) {
		startY = round(width * i/kafons)
		startX = 0
		
		xPos = startX
		yPos = startY
		
		platesLeft = plates
		
		for(j in 1:depth) {
			if(platesLeft == 0) {
				platedOut = platedOut + 1
				break;	
			}
			
			# 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) {
					dotCol = "gray"
				} else {
					covered[rndY,xPos] = 1
					dotCol = "blue"
				}
			
				# Did we hit a mine? If so make it red
				if(mineMatrix[rndY,xPos] == 1) {
					if(plotKafons) {
						points(xPos, rndY, col="red", pch=20, cex=0.5)
					}
					
					mineMatrix[rndY,xPos] = 0 # Mine explodes
					platesLeft = platesLeft - 1
				} else {
					if(plotKafons) {
						points(xPos, rndY, col=dotCol, pch=".")
					}
				}
			}
		}
	}
	
	history.numbKafons[trial] = k
	history.numbMines[trial] = startMines
	history.minesLeft[trial] = sum(mineMatrix)
	history.totalAreaCovered[trial] = sum(covered)
	history.platedOut[trial] = platedOut
	
	trial = trial + 1
}

# Plots of the results
# Mines left vs kafons run
# plot(history.numbKafons, (100*history.minesLeft/history.numbMines), xlab="Number of Kafons run", ylab="Percent of mines left", pch=20, col="blue")

# Area covered vs kafons run
# plot(history.numbKafons, 100*history.totalAreaCovered/(width*depth), xlab="Number of Kafons run", ylab="Percent of total area covered", pch=20, col="blue")

# Number of overlaps vs kafons run
# plot(history.numbKafons, (depth*history.numbKafons - history.totalAreaCovered) , xlab="Number of Kafons run", ylab="Amount of overlapped territory", pch=20, col="blue")

# Cost to destroy mines
# plot(history.numbKafons, ( (costPerKafon * history.numbKafons)/(history.numbMines - history.minesLeft) ) , xlab="Number of Kafons run", ylab="Cost per destroyed mine", pch=20, col="blue")