# This file should be named Data.R
# This script contains the data used in the assessment.
# For some of the parameters, the script composes arrays, in which 
# dimension 1 = years, dimension 2 = regions, dimension 3 = iterations.

##########
# The number of Monochamus adults sampled per inspection site
n_Monochamus_data <- function(n_y,n_i){
  A <- array(0,dim=c(n_y,15,n_i))
  a <- matrix(0,n_y,15)
  a[13,] = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
  a[14,] = c(0,0,0,0,0,0.83,0,0,0,0,0,0,0,0,0)
  a[15,] = c(0,0,0,0,0,1.5,0,0,0,0,0,0,0,9,0)
  a[16,] = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
  a[17,] = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
  a[18,] = c(0.5,0,0,0,0,6,2,1,0,0,0,0,0,0,0)
  a[19,] = c(0,0,0,0,0,1,0,0,0,0,0,0,4,0,0)
  for(i in 1:n_y){
    for(j in 1:15){
      A[i,j,] = a[i,j]
    }
  }
  return(A)} 

##########
# The number of inspection sites in the wood sampling component of the survey
N_wood_data <- function(n_y,n_i){
  A <- array(0,dim=c(n_y,15,n_i))
  a <- matrix(0,n_y,15)
  a[1,] = c(100,75,75,50,50,100,50,50,75,50,50,75,75,50,50)
  a[2,] = c(4,3,0,0,0,8,1,0,0,0,0,0,0,0,0)
  a[3,] = c(5,3,0,0,0,0,0,0,0,0,6,0,0,0,0)
  a[4,] = c(37,30,50,0,0,72,0,0,50,0,0,0,51,0,4)
  a[5,] = c(57,50,50,0,0,98,0,0,50,0,0,48,54,0,0)
  a[6,] = c(61,50,50,0,0,95,0,0,38,16,0,0,51,0,0)
  a[7,] = c(65,51,40,0,0,27,0,0,52,52,0,0,50,16,0)
  a[8,] = c(55,29,46,25,25,23,54,24,55,52,25,50,54,28,0)
  a[9,] = c(52,45,50,25,25,33,52,23,23,53,25,50,50,15,13)
  a[10,] = c(45,50,50,25,25,40,53,26,38,51,27,50,53,27,23)
  a[11,] = c(45,45,42,19,20,36,46,28,67,45,21,45,45,20,23)
  a[12,] = c(45,57,43,18,20,40,46,20,33,47,20,45,45,20,5)
  a[13,] = c(49,45,47,22,19,22,45,17,29,29,20,45,45,19,15)
  a[14,] = c(47,45,43,19,21,29,45,10,39,9,20,44,44,28,20)
  a[15,] = c(45,46,45,20,20,41,45,20,47,20,20,45,45,20,16)
  a[16,] = c(20,46,45,20,5,22,45,10,16,20,20,45,45,9,20)
  a[17,] = c(24,46,45,20,20,30,45,20,22,15,20,45,20,9,20)
  a[18,] = c(30,45,42,22,20,31,50,7,19,20,20,45,17,9,20)
  a[19,] = c(13,20,45,21,20,33,23,29,22,20,20,45,16,5,20)
  for(i in 1:n_y){
    for(j in 1:15){
      A[i,j,] = a[i,j]
    }
  }
  return(A)}

##########
# The number of inspection sites in the Monochamus trapping component of the survey
N_Monochamus_data <- function(n_y,n_i){
  A <- array(0,dim=c(n_y,15,n_i))
  a <- matrix(0,n_y,15)
  a[13,] = c(0,0,0,0,0,6,0,0,2,0,0,0,0,2,0)
  a[14,] = c(0,0,0,0,0,6,0,0,2,0,0,0,0,2,0)
  a[15,] = c(0,0,0,0,0,6,0,0,2,0,0,0,0,2,0)
  a[16,] = c(0,0,0,0,0,6,0,0,2,0,0,0,0,2,0)
  a[17,] = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
  a[18,] = c(2,0,0,0,0,6,2,1,0,0,0,0,0,0,0)
  a[19,] = c(2,0,0,0,0,1,0,0,0,0,0,0,4,0,0)
  for(i in 1:n_y){
    for(j in 1:15){
      A[i,j,] = a[i,j]
    }
  }
  return(A)}

##########
# The number of dead wood objects suitable for sampling per inspection site
p_wood_data <- function(n_y,n_i){
  
  # The number of Monochamus-suitable dead wood objects per km2 
  obj <- round(rpert(n_i,166,288,398,1))
  
  # The proportion of Monochamus-suitable dead wood objects that is suitable for sampling 
  psam <- runif(n_i,0.05,0.95)
  
  # The density of wood objects suitable for sampling, number/km2 
  D_wood = round(obj*psam)
  
  # The number of wood objects suitable for sampling per inspection site  
  p_wood_a = round(D_wood*site_wood)
  p_wood_a[p_wood_a < ceiling(1/DP_wood)] <- ceiling(1/DP_wood)
  
  A <-array(0,dim=c(n_y,15,n_i))
  for(i in 1:n_y){
    for(j in 1:15){
      A[i,j,] = p_wood_a
    }
  }
  return(A)}

##########
# The number of Monochamus adults per inspection site 
p_Monochamus_data <- function(n_y,n_i){
  
  # The number of dead wood objects occupied by Monochamus per km2
  obju <- round(rpert(n_i,13.28,28.8,47.76,1))
  
  # The number of Monochamus eggs laid per Monochamus-suitable dead-wood object 
  fobj <- rpert(n_i,6,31,88,1)
  
  # The proportion of Monochamus surviving from egg to egg-laying adults
  surv <- rpert(n_i,0.1,0.25,0.4,1)
  
  # The density of Monochamus adults, number/km2
  D_Monochamus = obju*fobj*surv
  
  # The number of Monochamus individuals per inspection site 
  p_Monochamus_a = round(D_Monochamus*site_Monochamus)
  p_Monochamus_a[p_Monochamus_a < ceiling(1/DP_Monochamus)] <- ceiling(1/DP_Monochamus)
  
  A <-array(0,dim=c(n_y,15,n_i))
  for(i in 1:n_y){
    for(j in 1:15){
      A[i,j,] = p_Monochamus_a
    }
  }
  return(A)}

########## 
# The area of entry sites (i.e. harbors, industrial areas and landfills) in the regions, km2
EA_data <- c(940,604,474,486,546,604,214,454,293,384,380,443,667,209,538)

##########
# The area of the target population, i.e. risk areas (i.e. areas with PWN host plants at 5 km radius from entry sites) in the regions, km2
Pop_data <- c(5260,5678,5486,7199,9475,7774,11102,12657,12301,13681,9423,8437,22717,13627,31229)

##########
# The mean area covered with PWN host plants in hypothetical PWN infestations with a 20-km diameter in the regions, km2
InfA_data <- c(156.0996,149.3567,185.0010,192.3717,201.9278,169.6825,201.5093,208.9255,200.0608,224.4638,210.4504,187.2985,211.6547,239.2596,209.7427)

##########
# The relative probability of invasion in the regions
RP_data <- function(n_y,n_i){
  A <- array(0,dim=c(n_y,15,n_i))
  for(i in 1:n_y){
    for(j in 1:n_i){
      A[i,,j] = EA_data/sum(EA_data)
    }
  }
  return(A)}

##########
# The effective probability of infestation for the import-export survey
# and the regional-level design prevalence for the early detection survey
DPr_adj_data <- function(scenario,n_y,n_i){
  
  # The relative probability of invasion in the regions
  RP = EA_data/sum(EA_data)
  
  # The proportion of the target population in the regions
  PropPop = Pop_data/sum(Pop_data)
  
  # Weighted probability of invasion in the regions
  WR = RP/sum(PropPop*RP)
  
  # Effective probability of infestation for the import-export survey
  A <- array(0,dim=c(n_y,15,n_i))
  if(scenario == 1){
    for(i in 1:n_y){
      for(j in 1:n_i){
        A[i,,j] = DPr*WR
      }
    }
  # Regional-level design prevalence for the early detection survey
  }else{
    RP_array <- array(0,dim=c(n_y,15,n_i))
    for(i in 1:n_y){
      for(j in 1:n_i){
        A[i,,j] = InfA_data/Pop_data
      }
    }
  }
  return(A)}

##########