# This file should be named Probability_of_freedom_function.R
# A function that returns the probability of freedom achieved by 2018
# for a range of probabilities of invasion for all the regions and Finland in one array 
# dimension 1 = years, dimension 2 = regions (regions columns 1-15, Finland column 16), dimension 3 = iterations

Probability_of_freedom_2018 <- function(scenario, DP_wood, DP_Monochamus, DPr, site_wood, site_Monochamus, Pinv_FI, n_i, n_y, PrioPfree_1){
  
  ##########
  # CALLING FOR THE DATA NEEDED
  
  # The number of wood objects suitable for sampling per inspection site
  p_wood = p_wood_data(n_y,n_i)
  
  # The number of Monochamus adults per inspection site
  p_Monochamus = p_Monochamus_data(n_y,n_i)
  
  # The number of wood samples collected per inspection site
  n_wood <- round(rpert(n_i,1,2,5,1))
  
  # The number of inspected sites in the wood sampling component of the survey
  N_wood = N_wood_data(n_y,n_i)
  
  # The number of Monochamus samples collected per inspection site
  n_Monochamus = n_Monochamus_data(n_y,n_i)
  
  # The number of inspected sites in the Monochamus trapping component of the survey
  N_Monochamus = N_Monochamus_data(n_y,n_i)
  
  # The relative probability of invasion in the regions
  RP = RP_data(n_y,n_i)
  
  # The effective probability of invasion for the import-export survey and 
  # the regional-level design prevalence for the early detection survey 
  DPr_adj = DPr_adj_data(scenario,n_y,n_i)
  
  ##########
  # CALCULATING INSPECTION SENSITIVITY
  
  ISe_wood = 1 - (1 - ((n_wood)/(p_wood - 0.5*(p_wood*DP_wood-1))))^(p_wood*DP_wood)
  ISe_Monochamus = 1 - (1 - ((n_Monochamus)/(p_Monochamus - 0.5*(p_Monochamus*DP_Monochamus-1))))^(p_Monochamus*DP_Monochamus)
  
  ##########
  # CALCULATING THE SENSITIVITY OF ANNUAL SURVEYS IN THE REGIONS
  # (separately for wood sampling and Monochamus trapping)
  
  GSe_w  = 1 - (1 - (DPr_adj*ISe_wood))^N_wood
  GSe_Monochamus = 1 - (1 - (DPr_adj*ISe_Monochamus))^N_Monochamus
  
  #########
  # CALCULATING THE SENSITIVITY OF ANNUAL SURVEYS (wood sampling and Monochamus trapping combined) 
  # AND THE PROBABILITY OF FREEDOM ACHIEVED BY THE DIFFERENT YEARS (2001-2018)
  
  # Arrays for storing the results 
  # w = wood, m = Monochamus, wm = wood and Monochamus
  GSe_wm <-array(0,dim=c(n_y,15,n_i))
  SSe_w_FI <-array(0,dim=c(n_y,n_i))
  SSe_m_FI <-array(0,dim=c(n_y,n_i))
  SSe_wm_FI <-array(0,dim=c(n_y,n_i))
  Pfree_w <-array(0,dim=c(n_y,15,n_i))
  Pfree_wm <-array(0,dim=c(n_y,15,n_i))
  Pfree_adj <-array(0,dim=c(n_y,15,n_i))
  Pfree_w_FI <-array(0,dim=c(n_y,n_i))
  Pfree_wm_FI <-array(0,dim=c(n_y,n_i))
  Pfree_adj_FI <-array(0,dim=c(n_y,n_i))
  Pfree_all <- array(0,dim=c(length(Pinv_FI),16,n_i))
  
  for (k in 1:n_i){
    for(t in 1:n_y){
      
      # Regions
      GSe_wm[t,,k] = 1 - (1-GSe_w[t,,k])*(1-GSe_Monochamus[t,,k])
      
      # Finland
      # Import-export
      if(scenario == 1){
        SSe_w_FI[t,k] = 1 - prod(1-GSe_w[t,,k])
        SSe_m_FI[t,k] = 1 - prod(1-GSe_Monochamus[t,,k])
        SSe_wm_FI[t,k] = 1 - (1-SSe_w_FI[t,k])*(1-SSe_m_FI[t,k])
        # Early detection
      }else{
        SSe_wm_FI[t,k] = sum(GSe_wm[t,,k]*RP[t,,k])
      }
    }
  }
  
  for (h in 1:length(Pinv_FI)){
    
    Pinv_country = Pinv_FI[h]
    
    # The probability of invasion in the regions
    Pinv_regions = RP*Pinv_country
    
    for(t in 1:n_y){
      
      if(t==1){  
        
        # Probability of freedom, Finland
        Pfree_wm_FI[t,] = PrioPfree_1 / (PrioPfree_1 + ((1-PrioPfree_1)*(1-SSe_wm_FI[t,])))
        
        # Adjusted probability of freedom, Finland
        Pfree_adj_FI[t,] = 1 - ( (1-Pfree_wm_FI[t,]) + Pinv_country - ((1-Pfree_wm_FI[t,])*Pinv_country) )
        
        # Probability of freedom, Regions
        Pfree_w[t,,] = PrioPfree_1 / (PrioPfree_1 + ((1-PrioPfree_1)*(1-GSe_w[t,,])))
        Pfree_wm[t,,] = Pfree_w[t,,] / (Pfree_w[t,,] + ((1-Pfree_w[t,,])*(1-GSe_Monochamus[t,,])))
        
        # Adjusted probability of freedom, Regions
        Pfree_adj[t,,] = 1 - ( (1-Pfree_wm[t,,]) + Pinv_regions[t,,] - ((1-Pfree_wm[t,,])*Pinv_regions[t,,]) )
        
      }else{
        
        # Probability of freedom, Finland
        Pfree_wm_FI[t,] = Pfree_adj_FI[t-1,] / (Pfree_adj_FI[t-1,] + ((1-Pfree_adj_FI[t-1,])*(1-SSe_wm_FI[t,])))
        
        # Adjusted probability of freedom, Finland
        Pfree_adj_FI[t,] = 1 - ( (1-Pfree_wm_FI[t,]) + Pinv_country - ((1-Pfree_wm_FI[t,])*Pinv_country) )
        
        # Probability of freedom, Regions
        Pfree_w[t,,] = Pfree_adj[t-1,,] / (Pfree_adj[t-1,,] + ((1-Pfree_adj[t-1,,])*(1-GSe_w [t,,])))
        Pfree_wm[t,,] = Pfree_w[t,,] / (Pfree_w[t,,] + ((1-Pfree_w[t,,])*(1-GSe_Monochamus[t,,])))
        
        # Adjusted probability of freedom, Regions
        Pfree_adj[t,,] = 1 - ( (1-Pfree_wm[t,,]) + Pinv_regions[t,,] - ((1-Pfree_wm[t,,])*Pinv_regions[t,,]) )
        
      }
    }
    
    ##########
    # RETURN ALL THE RESULTS IN ONE ARRAY
    
    # Regions in columns 1-15, Finland in column 16
    Pfree_all[h,1:15,] =  Pfree_wm[n_y,,]
    Pfree_all[h,16,] = Pfree_wm_FI[n_y,]
    
  }
  
  return(Pfree_all)}