Chessboard Model

library(ggplot2);library(magrittr); library(ggsci); library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(reshape2)

Rules (From McCue, 2005)

Ocean is an 8 x 5 “chessboard.” • Squares are 250 nm on a side. • A 12-boat wolf pack can search a square in a day. • The steps: 1. Update the wolf packs’ score. 2. Roll a die for each ongoing sighting, using theappropriate “rule of13rd” described below, to see ifit is converted into an attack. 3. Roll a die each convoy that is under attack to see if it gets rescued and returns to being unsighted, usingthe appropriate “rule of1/3rd” described below. 4. The SIGINT step: Assign a wolf pack to convergeupon any compromised convoy that is not yetthreatened. Then roll the die for ULTRA, if applica-ble, using the appropriate “rule of13rd” describedbelow, to return the convoy to unsighted status andthe corresponding pack to patrolling status. 5. Move the convoys as described below. 6. Move the packs as described below. 7. Check for new sightings. Any unsighted or threat-ened convoy in the same square as a wolf packbecomes sighted, including ones that just becameunsighted as a result of ULTRA

ConvoyTool = function(Width = 5, Depth = 8, phold = .66, DaysToMove = 2){

ConvoyMat = matrix(0, nrow = Width, ncol = Depth)
SubMat = matrix(0, nrow = Width, ncol = Depth)
Score = 0
turn = 0

#initialize Positions

StartX = 1
StartY = 1
ConvoyMat[StartX, StartY] = 1

SubX = sample(1:Width, 1)
SubY = sample(1:Depth, 1)

SubMat[SubX, SubY] = 1

#Score adjudicated when Wolfpack attacks a convoy
ConD = 0

### Loop over the game

while(ConD <= Depth){
ConvoyOld = ConvoyMat #Save a copy for below

if(which(ConvoyMat == 1) == which(SubMat == 1)){Score %<>% `+`(1)}

turn %<>% `+`(1)
ConW = which(ConvoyMat == 1, arr.ind = TRUE)[1]
ConD = which(ConvoyMat == 1, arr.ind = TRUE)[2]

#Move Convoy on Even Turns only 
if(turn %% DaysToMove == 0){
  ConD = ConD + 1
  
  ConW = if(ConW == 1){ConW = 2
  }else if(ConW == Width){
    ConW = Width - 1
    }else if(runif(1) < .5){
      ConW = ConW - 1
      }else{
        ConW = ConW + 1
      }
}



ConvoyMat[,] = 0

ConvoyMat[ConW, min(ConD, Depth)] = 1 #Min Statemetn prevents crashes

#Move Hunters

SubW = 1
SubD = 1 #Error generation for base cases

SubWOld = which(SubMat == 1, arr.ind = TRUE)[1]
SubDOld = which(SubMat == 1, arr.ind = TRUE)[2]

#If in contact, 2/3 chance to stay in contact
if(which(ConvoyOld == 1) == which(SubMat == 1)){
  if(runif(1) < phold){
    SubW = ConW
    SubD = ConD
  }
}else{ #Move the Hunters
 if(SubWOld == 1){SubW = 2} else if(SubWOld == Width){
    SubW = Width - 1} else if(runif(1) < .5){SubW = SubWOld -1}else{
      SubW = SubWOld + 1
    }

  if(SubDOld == 1){SubD =2}else if(SubDOld == Depth){
    SubD = Depth - 1}else if(runif(1) < .5){SubD = SubDOld -1} else{
      SubD = SubDOld + 1
    }
}

SubMat[,] = 0
SubMat[SubW, min(Depth, SubD)] = 1
}
return(Score)
}

ConvoyMaster = function(Nreps = 100, phold = .66, DaysToMove = 2){
  Vec = rep(-1, length = Nreps)
  for(i in 1:Nreps){
    Vec[i] = ConvoyTool(Width = 5, Depth = 8, phold = phold, DaysToMove = DaysToMove)
  }
  return(Vec)
}

set.seed(2255)
Basic = ConvoyMaster(Nreps = 100) #Actual production had 10000 reps

sum(Basic==0)

## [1] 72

Contact = Basic[Basic>0]

theme_set(theme_minimal())

data.frame(Contact = Contact) %>% ggplot(aes(x = Contact)) + geom_density(color = 'darkblue') + ggtitle("Contact Days") + labs(subtitle = "Basic Convoy Model") + xlab("Contact Days") + scale_y_continuous(labels = scales::percent)

High Speed Excursion

set.seed(2255)
HighSpeed = ConvoyMaster(Nreps = 100, phold = .8, DaysToMove = 1)

#Let’s try starting the submarine closer to the convoy

ConvoyTool2 = function(Width = 5, Depth = 8, phold = .66, DaysToMove = 2){

ConvoyMat = matrix(0, nrow = Width, ncol = Depth)
SubMat = matrix(0, nrow = Width, ncol = Depth)
Score = 0
turn = 0

#initialize Positions

StartX = 1
StartY = 1
ConvoyMat[StartX, StartY] = 1

SubX = sample(1:3, 1)
SubY = sample(1:3, 1)
SubMat[SubX, SubY] = 1

#Score adjudicated when Wolfpack attacks a convoy
ConD = 0

### Loop over the game

while(ConD <= Depth){
ConvoyOld = ConvoyMat #Save a copy for below

if(which(ConvoyMat == 1) == which(SubMat == 1)){Score %<>% `+`(1)}

turn %<>% `+`(1)
ConW = which(ConvoyMat == 1, arr.ind = TRUE)[1]
ConD = which(ConvoyMat == 1, arr.ind = TRUE)[2]

#Move Convoy on Even Turns only 
if(turn %% DaysToMove == 0){
  ConD = ConD + 1
  
  ConW = if(ConW == 1){ConW = 2
  }else if(ConW == Width){
    ConW = Width - 1
    }else if(runif(1) < .5){
      ConW = ConW - 1
      }else{
        ConW = ConW + 1
      }
}



ConvoyMat[,] = 0

ConvoyMat[ConW, min(ConD, Depth)] = 1 #Min Statemetn prevents crashes

#Move Hunters

SubW = 1
SubD = 1 #Error generation for base cases

SubWOld = which(SubMat == 1, arr.ind = TRUE)[1]
SubDOld = which(SubMat == 1, arr.ind = TRUE)[2]

#If in contact, 2/3 chance to stay in contact
if(which(ConvoyOld == 1) == which(SubMat == 1)){
  if(runif(1) < phold){
    SubW = ConW
    SubD = ConD
  }
}else{ #Move the Hunters
 if(SubWOld == 1){SubW = 2} else if(SubWOld == Width){
    SubW = Width - 1} else if(runif(1) < .5){SubW = SubWOld -1}else{
      SubW = SubWOld + 1
    }

  if(SubDOld == 1){SubD =2}else if(SubDOld == Depth){
    SubD = Depth - 1}else if(runif(1) < .5){SubD = SubDOld -1} else{
      SubD = SubDOld + 1
    }
}

SubMat[,] = 0
SubMat[SubW, min(Depth, SubD)] = 1
}
return(Score)
}

ConvoyMaster2 = function(Nreps = 100, phold = .66, DaysToMove = 2){
  Vec = rep(-1, length = Nreps)
  for(i in 1:Nreps){
    Vec[i] = ConvoyTool2(Width = 5, Depth = 8, phold = phold, DaysToMove = DaysToMove)
  }
  return(Vec)
}

set.seed(2255)

CloseStart = ConvoyMaster2(Nreps = 100, phold = .5, DaysToMove = 2)