From 1f38c1a3e7863200006a3d7c366a3397de947aac Mon Sep 17 00:00:00 2001
From: David Pleydell <david.pleydell@inrae.fr>
Date: Wed, 19 May 2021 23:04:51 +0200
Subject: [PATCH 1/3] Started verbose_plotting.R

---
 nimble/model12/verbose_plotting.R | 38 +++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 nimble/model12/verbose_plotting.R

diff --git a/nimble/model12/verbose_plotting.R b/nimble/model12/verbose_plotting.R
new file mode 100644
index 0000000..36a49ab
--- /dev/null
+++ b/nimble/model12/verbose_plotting.R
@@ -0,0 +1,38 @@
+library(here)
+SCENARIO <- 1 # 4
+(modelDir = here::here("nimble/model12"))         # For selecting model version
+(mcmcDir  = paste0(modelDir,"/MCMCall1000"))      # For selecting samples.csv. Can be a symbolic link to a different directory if simualtion code is more recent that estimation code (e.g. fence added since MCMC).
+(simDir   = paste0(mcmcDir,"/simsCases_Verbose")) # For selecting model version
+
+setwd(simDir)
+simFiles = sort(dir()[grep(paste0("scenario",SCENARIO), dir())])
+
+tStop = 230
+nSim  = length(simFiles)
+
+nSim = 100
+## E, I, C globally
+par(mfrow=n2mfrow(3))
+tRange = -30:tStop
+plot(tRange, tRange, typ="n", ylab="E", xlab="Day", ylim=c(0,2200))
+for (sim in 1:nSim) { # sim=1
+  if (sim %% round(nSim/10) == 1) print(sim)
+  load(simFiles[sim])
+  lines(tVec, rowSums(E), col=rgb(0,0,0,0.1))
+}
+plot(tRange, tRange, typ="n", ylab="I", ylim=c(0,1800))
+for (sim in 1:nSim) { # sim=1
+  if (sim %% round(nSim/10) == 1) print(sim)
+  load(simFiles[sim])
+  lines(tVec, rowSums(I), col=rgb(0,0,0,0.1))
+}
+plot(tRange, tRange, typ="n", ylab="C", ylim=c(0,10000))
+for (sim in 1:nSim) { # sim=1
+  if (sim %% round(nSim/10) == 1) print(sim)
+  load(simFiles[sim])
+  lines(tVec, rowSums(C), col=rgb(0,0,0,0.1))
+}
+
+
+## Can the output be sparse? i.e. avoid pixels that do nothing?
+# plot(colSums(outputList[[sim]]$states[5,,]))
-- 
GitLab


From c5acee258795ea866bfdf2e4e2861ef6fd089176 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Thu, 20 May 2021 11:37:57 +0200
Subject: [PATCH 2/3] Edits in verbose_plotting.R

---
 nimble/model12/verbose_plotting.R             | 162 +++++++++++++++---
 .../model12/verbose_simWildBoarFromMCMC12.R   |   4 +-
 2 files changed, 138 insertions(+), 28 deletions(-)

diff --git a/nimble/model12/verbose_plotting.R b/nimble/model12/verbose_plotting.R
index 36a49ab..81714c6 100644
--- a/nimble/model12/verbose_plotting.R
+++ b/nimble/model12/verbose_plotting.R
@@ -1,38 +1,148 @@
-library(here)
-SCENARIO <- 1 # 4
-(modelDir = here::here("nimble/model12"))         # For selecting model version
-(mcmcDir  = paste0(modelDir,"/MCMCall1000"))      # For selecting samples.csv. Can be a symbolic link to a different directory if simualtion code is more recent that estimation code (e.g. fence added since MCMC).
-(simDir   = paste0(mcmcDir,"/simsCases_Verbose")) # For selecting model version
+## rm(list=ls())
+library(sf)
+library(tmap)
+library(tmaptools)
+library(dplyr)
 
-setwd(simDir)
+
+## library(here)
+## (modelDir = here::here("nimble/model12"))         # For selecting model version
+## (mcmcDir  = paste0(modelDir,"/MCMCall1000"))      # For selecting samples.csv. Can be a symbolic link to a different directory if simualtion code is more recent that estimation code (e.g. fence added since MCMC).
+## (simDir   = paste0(mcmcDir,"/simsCases_Verbose")) # For selecting model version
+## setwd(simDir)
+
+## Choose simulation files for a given SCENARIO
+SCENARIO <- 1 # could also be 2 (and eventually 3, 4, 5, 6, 7 or 8)
 simFiles = sort(dir()[grep(paste0("scenario",SCENARIO), dir())])
 
+## Set two constants
 tStop = 230
 nSim  = length(simFiles)
 
-nSim = 100
-## E, I, C globally
-par(mfrow=n2mfrow(3))
+## Plot trajectories of 3 hidden states and 3 observed states (detected cases)
+nSim4plot = 10
 tRange = -30:tStop
-plot(tRange, tRange, typ="n", ylab="E", xlab="Day", ylim=c(0,2200))
-for (sim in 1:nSim) { # sim=1
-  if (sim %% round(nSim/10) == 1) print(sim)
-  load(simFiles[sim])
-  lines(tVec, rowSums(E), col=rgb(0,0,0,0.1))
+ALPHA = 0.2
+par(mfrow=n2mfrow(6))
+## Exposed
+plot(tRange, tRange, typ="n", ylab="Exposed", xlab="Day", ylim=c(0,2200))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(E), col=rgb(0,0,0,ALPHA))
+}
+## Infectious
+plot(tRange, tRange, typ="n", ylab="Infectious", xlab="Day", ylim=c(0,1800))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(I), col=rgb(0,0,0,ALPHA))
+}
+## Carcass
+plot(tRange, tRange, typ="n", ylab="Carcass", xlab="Day", ylim=c(0,10000))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(C), col=rgb(0,0,0,ALPHA))
 }
-plot(tRange, tRange, typ="n", ylab="I", ylim=c(0,1800))
-for (sim in 1:nSim) { # sim=1
-  if (sim %% round(nSim/10) == 1) print(sim)
-  load(simFiles[sim])
-  lines(tVec, rowSums(I), col=rgb(0,0,0,0.1))
+## Hunted & tested +ve
+plot(tRange, tRange, typ="n", ylab="Hunted & Tested +ve", xlab="Day", ylim=c(0,200))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(HP), col=rgb(0,0,0,ALPHA))
 }
-plot(tRange, tRange, typ="n", ylab="C", ylim=c(0,10000))
-for (sim in 1:nSim) { # sim=1
-  if (sim %% round(nSim/10) == 1) print(sim)
-  load(simFiles[sim])
-  lines(tVec, rowSums(C), col=rgb(0,0,0,0.1))
+## Active search
+plot(tRange, tRange, typ="n", ylab="Active search", xlab="Day", ylim=c(0,50))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(AS), col=rgb(0,0,0,ALPHA))
 }
+## Passive search
+plot(tRange, tRange, typ="n", ylab="Passive search", xlab="Day", ylim=c(0,25))
+abline(h=0)
+for (sim in 1:nSim4plot) { # sim=1
+  if (sim %% round(nSim4plot/10) == 1) print(sim)
+  load(simFiles[sample(nSim, 1)]) # Choose one file at random & load it
+  lines(tVec, rowSums(PS), col=rgb(0,0,0,ALPHA))
+}
+
 
+#############
+## Mapping ##
+#############
+#drake::loadd(hexAll, hexSub_at_D110, fence, huntZone, admin)
+#save(hexAll, hexSub_at_D110, fence, huntZone, admin, file="objectsForMapping.Rdata")
+load("objectsForMapping.Rdata")
 
-## Can the output be sparse? i.e. avoid pixels that do nothing?
-# plot(colSums(outputList[[sim]]$states[5,,]))
+## Choose one file at random & load it
+load(simFiles[sample(nSim, 1)])
+## Select simulated data from one given day
+DAY = 80
+ZOOM = TRUE
+hexAll$E  = E[which(tVec==DAY),]
+hexAll$I  = I[which(tVec==DAY),]
+hexAll$C  = C[which(tVec==DAY),]
+hexAll$HP = HP[which(tVec==DAY),]
+hexAll$AS = AS[which(tVec==DAY),]
+hexAll$PS = PS[which(tVec==DAY),]
+## Define map extent
+if (ZOOM==TRUE) {
+  mapBase = tm_shape(hexSub_at_D110) + tm_polygons(col = "pForest", lwd=0.2, border.col="grey", palette=get_brewer_pal("Greens", n = 10, plot=FALSE))
+} else {
+  mapBase = tm_shape(hexAll) + tm_polygons(col = "pForest", lwd=0.2, border.col="grey", palette=get_brewer_pal("Greens", n = 10, plot=FALSE))
+}
+## Define finishing overlays
+mapOverays =
+  tm_shape(admin) + tm_borders(lwd=2) +
+  tm_shape(fence) + tm_borders(col="darkviolet", lwd=2) +
+  tm_shape(huntZone) + tm_borders(col="red", lwd=2)
+## Map of exposed
+mapE = mapBase
+if(any(hexAll$E>0)) {
+  mapE = mapE + tm_shape(hexAll %>% filter(hexAll$E > 0)) + tm_polygons(col="E", border.alpha=1)
+}
+mapE = mapE + mapOverays
+## Map of infectious
+mapI = mapBase
+if(any(hexAll$I>0)) {
+  mapI = mapI +
+    tm_shape(hexAll %>% filter(hexAll$I > 0)) + tm_polygons(col="I", border.alpha=1)
+}
+mapI = mapI + mapOverays
+## Map of Carcasses
+mapC = mapBase
+if(any(hexAll$C>0)) {
+  mapC = mapC +
+    tm_shape(hexAll %>% filter(hexAll$C > 0)) + tm_polygons(col="C", border.alpha=1)
+}
+mapC = mapC + mapOverays
+## Map of Hunted & tested +ve
+mapHP = mapBase
+if(any(hexAll$HP>0)) {
+  mapHP = mapHP +
+    tm_shape(hexAll %>% filter(hexAll$HP > 0)) + tm_polygons(col="HP", border.alpha=1)
+}
+mapHP = mapHP + mapOverays
+## Map of Active search
+mapAS = mapBase
+if(any(hexAll$AS>0)) {
+  mapAS = mapAS +
+    tm_shape(hexAll %>% filter(hexAll$AS > 0)) + tm_polygons(col="AS", border.alpha=1)
+}
+mapAS = mapAS + mapOverays
+## Map of Passive search
+mapPS = mapBase
+if(any(hexAll$PS>0)) {
+  mapPS = mapPS +
+    tm_shape(hexAll %>% filter(hexAll$PS > 0)) + tm_polygons(col="PS", border.alpha=1)
+}
+mapPS = mapPS + mapOverays
+## Plot all maps
+tmap_arrange(mapE, mapI, mapC, mapHP, mapAS, mapPS)
diff --git a/nimble/model12/verbose_simWildBoarFromMCMC12.R b/nimble/model12/verbose_simWildBoarFromMCMC12.R
index 720452c..a925d6c 100644
--- a/nimble/model12/verbose_simWildBoarFromMCMC12.R
+++ b/nimble/model12/verbose_simWildBoarFromMCMC12.R
@@ -1,7 +1,7 @@
 # source(here::here("nimble/model12/verbose_simWildBoarFromMCMC12.R"))
 # rm(list=ls())
 
-SCENARIO <- 1 # 4
+SCENARIO <- 2 # 4
 
 # This script jumps through MCMC output and generates weighted means & stdevs of the simulations
 
@@ -166,7 +166,7 @@ cSimFromMCMC <- compileNimble(rSimFromMCMC)
 ##%%%%%%%%%%%%%%%%%%%%%
 ## Run simulations ####
 ##%%%%%%%%%%%%%%%%%%%%%
-nSim         = 1000 ## 10000 # 5000
+nSim         = 10000 ## 10000 # 5000
 tFence       = 60
 flagFence    = c(1, 0) ## if 0 -> set omegaFence to 0
 flagHuntZone = c(1, 0) ## if 0 -> set log_tauHArmy to -Inf (if it can handle it)
-- 
GitLab


From bcbd81fcc0295f45d7b0e899fd6de8bb580ced01 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Sun, 23 May 2021 23:38:19 +0200
Subject: [PATCH 3/3] Updates in verbose_plotting.R

---
 nimble/model12/verbose_plotting.R | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/nimble/model12/verbose_plotting.R b/nimble/model12/verbose_plotting.R
index 81714c6..a95f64d 100644
--- a/nimble/model12/verbose_plotting.R
+++ b/nimble/model12/verbose_plotting.R
@@ -4,16 +4,14 @@ library(tmap)
 library(tmaptools)
 library(dplyr)
 
-
-## library(here)
-## (modelDir = here::here("nimble/model12"))         # For selecting model version
-## (mcmcDir  = paste0(modelDir,"/MCMCall1000"))      # For selecting samples.csv. Can be a symbolic link to a different directory if simualtion code is more recent that estimation code (e.g. fence added since MCMC).
-## (simDir   = paste0(mcmcDir,"/simsCases_Verbose")) # For selecting model version
-## setwd(simDir)
+## Set the working directory to the directory storing the simulation output (Rdata) files
+setwd("/home/david/asf-challenge/nimble/model12/MCMCall1000/simsCases_Verbose")
 
 ## Choose simulation files for a given SCENARIO
-SCENARIO <- 1 # could also be 2 (and eventually 3, 4, 5, 6, 7 or 8)
+SCENARIO <- 3 # could also be 2 (and eventually 3, 4, 5, 6, 7 or 8)
 simFiles = sort(dir()[grep(paste0("scenario",SCENARIO), dir())])
+simFiles = simFiles[grep("Rdata", simFiles)]
+length(simFiles) ## Should be 10000
 
 ## Set two constants
 tStop = 230
-- 
GitLab