# 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")