# This file should be named Sensitivity.R
# Calculates the sensitivity of the annual surveys in 2000-2018 
# for all the regions and Finland and plots them as figures

##########
# 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("Sensitivity_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

##########

# 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

# Calling for the function that calculates the sensitivities
Scenario_1 = Sensitivity(1, DP_wood, DP_Monochamus, DPr, site_wood, site_Monochamus, n_i)

##########

# 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 IN THE CALCULATION BUT NEEDS TO BE GIVEN A VALUE 
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

# Calling for the function that calculates the sensitivities
Scenario_2 = Sensitivity(2, DP_wood, DP_Monochamus, DPr, site_wood, site_Monochamus, n_i)

##########
# FINDING THE 2.5, 50 AND 97.5 PERCENTILES

# Matrixes for storing the results
# Se1 = import-export (scenario 1)
# Se2 = early detection (scenario 2)
Se1_2.5 <- matrix(0,n_y,16)
Se1_50 <- matrix(0,n_y,16)
Se1_97.5 <- matrix(0,n_y,16)
Se2_2.5 <- matrix(0,n_y,16)
Se2_50 <- matrix(0,n_y,16)
Se2_97.5 <- matrix(0,n_y,16)

for(i in 1:n_y){
  for(j in 1:16){
    Se1_2.5[i,j] = quantile(Scenario_1[i,j,],0.025)
    Se1_50[i,j] = quantile(Scenario_1[i,j,],0.5)
    Se1_97.5[i,j] = quantile(Scenario_1[i,j,],0.975)
    Se2_2.5[i,j] = quantile(Scenario_2[i,j,],0.025) 
    Se2_50[i,j] = quantile(Scenario_2[i,j,],0.5) 
    Se2_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.4
TitleSize<-0.8 
AxisSize<-0.75
Xlabels<-c(2000,2006,2012,2018)
Ylabels<-c(0.00,0.20,0.40,0.60,0.80,1.00)
AScenCol<-rgb(255,0,0,max=255,alpha=100)
BScenCol<-rgb(0,0,255,max=255,alpha=100)
APolCol<-rgb(255,0,0,max=255,alpha=100)
BPolCol<-rgb(0,0,255,max=255,alpha=100)
Whiskers<-0.01

plot(Y, Se1_50[,1], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0.00,1.00), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,1], Y, Se1_97.5[,1], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,1], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,1], Y, Se2_97.5[,1], col=BScenCol,length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,2], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,2], Y,Se1_97.5[,2], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,2], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,2], Y,Se2_97.5[,2], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,3], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,3], Y, Se1_97.5[,3], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,3], pch=16, col = BScenCol)
arrows(Y, Se2_2.5[,3], Y, Se2_97.5[,3], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,4], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,4], Y, Se1_97.5[,4], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,4], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,4], Y, Se2_97.5[,4], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,5], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,5], Y, Se1_97.5[,5], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,5], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,5], Y, Se2_97.5[,5], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,6], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,6], Y, Se1_97.5[,6], col = AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,6], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,6], Y, Se2_97.5[,6], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,7], type = "p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,7], Y, Se1_97.5[,7], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,7], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,7], Y, Se2_97.5[,7], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,8], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,8], Y, Se1_97.5[,8], col = AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,8], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,8], Y, Se2_97.5[,8], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,9], type="p", pch=16,col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,9], Y, Se1_97.5[,9], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,9], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,9], Y, Se2_97.5[,9], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,10], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,10], Y, Se1_97.5[,10], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,10], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,10], Y, Se2_97.5[,10], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,11], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,11], Y, Se1_97.5[,11], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,11], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,11], Y, Se2_97.5[,11], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,12], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,12], Y, Se1_97.5[,12], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,12], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,12], Y, Se2_97.5[,12], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,13], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,13], Y, Se1_97.5[,13], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,13], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,13], Y, Se2_97.5[,13], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,14], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,14], Y, Se1_97.5[,14], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,14], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,14], Y, Se2_97.5[,14], col=BScenCol, length=Whiskers, angle=90, code=3)
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(Y, Se1_50[,15], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="", ylab="")
arrows(Y, Se1_2.5[,15], Y, Se1_97.5[,15], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,15], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,15], Y, Se2_97.5[,15], col=BScenCol, length=Whiskers, angle=90, code=3)
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('Year', side=1, outer=TRUE, line=1.2, cex=0.75)
mtext('Sensitivity', 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(Y, Se1_50[,16], type="p", pch=16, col=AScenCol, frame.plot=FALSE, ylim=c(0,1), axes=FALSE, xlab="Year", ylab="Sensitivity")
arrows(Y, Se1_2.5[,16], Y, Se1_97.5[,16], col=AScenCol, length=Whiskers, angle=90, code=3)
points(Y, Se2_50[,16], pch=16, col=BScenCol)
arrows(Y, Se2_2.5[,16], Y, Se2_97.5[,16], col=BScenCol, length=Whiskers, angle=90, code=3)
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 sensitivity of the annual surveys. The dots denote the medians', adj=0, side=1, outer=FALSE, line=5, cex=1)
mtext('and the bars the 95% confidence intervals of the assessment results.', adj=0, side=1, outer=FALSE, line=6, cex=1)
mtext('Red denotes the import-export surveys and blue the early detection surveys.', adj=0, side=1, outer=FALSE, line=7, cex=1)