# This file should be named Probability_of_freedom.R
# This script calculates the probability of freedom achieved by 2018
# for a range of probabilities of invasion for all the regions and Finland and plots them as figures.
# If the number of iterations is high, running the script can take several minutes.
##########

# Installing and loading the package mc2d version 0.1-18
installed = rownames(installed.packages())
if (! "mc2d" %in% installed) {
  install.packages("mc2d")
}
library(mc2d)

# Loading the data and functions 
source("Data.R")
source("Probability_of_freedom_function.R")

n_i = 10000 # The number of iterations
Y = seq(2000,2018,1) # The studied years
n_y = length(Y) # The number of years studied
PrioPfree_1 = 0.5 # The initial prior probability of freedom
F_inv_FI = seq(2,100,1) # The mean time between invasions to Finland
Pinv_FI = 1/F_inv_FI # The probability of invasion to Finland in one year

##########
# IMPORT-EXPORT SURVEY
DP_wood = 0.12 # Local-level design prevalence, wood
DP_Monochamus = 0.09 # Local-level design prevalence, Monochamus
DPr = 0.01 # National-level design prevalence
site_wood = 0.35 # Inspection site size in the wood sampling component, km2
site_Monochamus = 0.25 # Inspection site size in the Monochamus trapping component, km2

Scenario_1 = Probability_of_freedom_2018(1, DP_wood, DP_Monochamus, DPr, site_wood, site_Monochamus, Pinv_FI, n_i, n_y, PrioPfree_1)

##########
# EARLY DETECTION SURVEY
DP_wood = 0.06 # Local-level design prevalence
DP_Monochamus = 0.045 # Local-level design prevalence, Monochamus
DPr = 0 # National-level design prevalence - NOT USED
site_wood = 0.63 # Inspection site size in the wood sampling component, km2
site_Monochamus = 0.25 # Inspection site size in the Monochamus trapping component, km2

Scenario_2 = Probability_of_freedom_2018(2, DP_wood, DP_Monochamus, DPr, site_wood, site_Monochamus, Pinv_FI, n_i, n_y, PrioPfree_1)

##########
# FINDING THE 2.5 AND 97.5 PERCENTILES

# Matrixes for storing the results
# PF1 = import-export (scenario 1)
# PF2 = early detection (scenario 2)
PF1_2.5 <- matrix(0,length(Pinv_FI),16)
PF1_97.5 <- matrix(0,length(Pinv_FI),16)
PF2_2.5 <- matrix(0,length(Pinv_FI),16)
PF2_97.5 <- matrix(0,length(Pinv_FI),16)

for(i in 1:length(Pinv_FI)){
  for(j in 1:16){
    PF1_2.5[i,j] = quantile(Scenario_1[i,j,],0.025) 
    PF1_97.5[i,j] = quantile(Scenario_1[i,j,],0.975)
    PF2_2.5[i,j] = quantile(Scenario_2[i,j,],0.025)   
    PF2_97.5[i,j] = quantile(Scenario_2[i,j,],0.975)     
  }
}

##########
# FIGURES
# The first figure presents results for the 15 administrative regions, 
# while the second figure presents respective results for Finland.  

par(mfrow=c(5,3))
par(pin=c(2,2))
par(mar=c(1,1,2,2))
par(mgp=c(1,0.8,0)) 
par(oma=c(2.5,3,0,0)) 

SetLine<--0.2
TitleSize<-0.8
AxisSize<-0.75
Xlabels<-c(0,10,20,30,40,50)
Ylabels<-c(0.00,0.20,0.40,0.60,0.80,1.00)
AScenCol<-rgb(0,0,0,max=255,alpha=0)
BScenCol<-rgb(0,0,0,max=255,alpha=0)
APolCol<-rgb(255,0,0,max=255,alpha=100)
BPolCol<-rgb(0,0,255,max=255,alpha=100)

plot(F_inv_FI, PF1_2.5[,1], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0.00,1.00), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,1], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,1], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,1], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,1], rev(PF1_97.5[,1])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,1], rev(PF2_97.5[,1])), col=BPolCol, border=NA)
title("Uusimaa", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,2], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,2], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,2], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,2], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,2], rev(PF1_97.5[,2])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,2], rev(PF2_97.5[,2])), col=BPolCol, border=NA)
title("Varsinais-Suomi", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,3], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,3], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,3], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,3], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,3], rev(PF1_97.5[,3])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,3], rev(PF2_97.5[,3])), col=BPolCol, border=NA)
title("Satakunta", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,4], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,4], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,4], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,4], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,4], rev(PF1_97.5[,4])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,4], rev(PF2_97.5[,4])), col=BPolCol, border=NA)
title("Häme", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,5], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,5], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,5], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,5], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,5], rev(PF1_97.5[,5])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,5], rev(PF2_97.5[,5])), col=BPolCol, border=NA)
title("Pirkanmaa", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,6], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,6], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,6], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,6], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,6], rev(PF1_97.5[,6])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,6], rev(PF2_97.5[,6])), col=BPolCol, border=NA)
title("Kaakkois-Suomi", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,7], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,7], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,7], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,7], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,7], rev(PF1_97.5[,7])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,7], rev(PF2_97.5[,7])), col=BPolCol, border=NA)
title("Etelä-Savo", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,8], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,8], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,8], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,8], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,8], rev(PF1_97.5[,8])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,8], rev(PF2_97.5[,8])), col=BPolCol, border=NA)
title("Pohjois-Savo", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,9], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,9], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,9], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,9], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,9], rev(PF1_97.5[,9])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,9], rev(PF2_97.5[,9])), col=BPolCol, border=NA)
title("Pohjois-Karjala", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,10], type="l", col=AScenCol, frame.plot=FALSE, ylim = c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,10], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,10], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,10], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,10], rev(PF1_97.5[,10])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,10], rev(PF2_97.5[,10])), col=BPolCol, border=NA)
title("Keski-Suomi", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,11], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,11], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,11], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,11], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,11], rev(PF1_97.5[,11])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,11], rev(PF2_97.5[,11])), col=BPolCol, border=NA)
title("Etelä-Pohjanmaa", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,12], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,12], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,12], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,12], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,12], rev(PF1_97.5[,12])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,12], rev(PF2_97.5[,12])), col=BPolCol, border=NA)
title("Pohjanmaa", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,13], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,13], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,13], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,13], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,13], rev(PF1_97.5[,13])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,13], rev(PF2_97.5[,13])), col=BPolCol, border=NA)
title("Pohjois-Pohjanmaa", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,14], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,14], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,14], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,14], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,14], rev(PF1_97.5[,14])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,14], rev(PF2_97.5[,14])), col=BPolCol, border=NA)
title("Kainuu", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

plot(F_inv_FI, PF1_2.5[,15], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="", ylab="")
lines(F_inv_FI, PF1_97.5[,15], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,15], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,15], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,15], rev(PF1_97.5[,15])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,15], rev(PF2_97.5[,15])), col=BPolCol, border=NA)
title("Lappi", line=SetLine, cex.main=TitleSize)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=AxisSize)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=AxisSize)

mtext('The mean time between invasions to Finland, years', side=1, outer=TRUE, line=1.2, cex=0.75)
mtext('The probability of freedom', side=2, outer=TRUE, line=1.5, cex=0.75)

par(mfrow=c(1,1))
par(pin=c(2,2))
par(mar=c(10,4,4,4))
par(mgp=c(2.7,1,0)) 
par(oma=c(1,0,0,0)) 

plot(F_inv_FI, PF1_2.5[,16], type="l", col=AScenCol, frame.plot=FALSE, ylim=c(0,1), xlim=c(0,50), axes=FALSE, xlab="The mean time between invasions to Finland, years", ylab="The probability of freedom")
lines(F_inv_FI, PF1_97.5[,16], col=AScenCol)
lines(F_inv_FI, PF2_2.5[,16], col=BScenCol)
lines(F_inv_FI, PF2_97.5[,16], col=BScenCol)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF1_2.5[,16], rev(PF1_97.5[,16])), col=APolCol, border=NA)
polygon(c(F_inv_FI,rev(F_inv_FI)), c(PF2_2.5[,16], rev(PF2_97.5[,16])), col=BPolCol, border=NA)
title("Finland", line=0)
axis(side=1, at=Xlabels, labels=Xlabels, cex.axis=1)
axis(side=2, at=Ylabels, labels=Ylabels, las=1, cex.axis=1)
mtext('The probability of freedom from PWN achieved by 2018. The colored areas', adj=0, side=1, outer=FALSE, line=6, cex=1)
mtext('show the 95% confidence intervals of the assessment results. Red denotes', adj=0, side=1, outer=FALSE, line=7, cex=1)
mtext('the import-export surveys, and blue denotes the early detection surveys.', adj=0, side=1, outer=FALSE, line=8, cex=1)