From 0c1f84dc37a3fb509022f7692496e9531547dd5a Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 04:58:04 +0100
Subject: [PATCH 1/9] Modifications for improved active search... but... we
don't have time to fit it. Zut!
---
nimble/model9/fitWildBoarModel9.R | 39 +++++++++++++++++-----------
nimble/model9/model_definition.R | 37 +++++++++++++++++++-------
nimble/model9/simWildBoarFromMCMC9.R | 39 ++++++++++++++++++----------
src/functions.R | 12 ++++++---
src/plan.R | 11 ++++++--
5 files changed, 94 insertions(+), 44 deletions(-)
diff --git a/nimble/model9/fitWildBoarModel9.R b/nimble/model9/fitWildBoarModel9.R
index 7526683..dcf581d 100644
--- a/nimble/model9/fitWildBoarModel9.R
+++ b/nimble/model9/fitWildBoarModel9.R
@@ -40,8 +40,12 @@ loadd(admin,
hexAll,
hexCentroidsSub,
fence,
- probNeighAS1km,
- probNeighAS2km,
+ weightTauAHome1km,
+ weightTauAAway1km,
+ weightTauAHome2km,
+ weightTauAAway2km,
+ ## probNeighAS1km,
+ ## probNeighAS2km,
outbreaks,
stepsPerChunkSub,
wildBoarArraySub, # observations aggregated over ten day intervals, hexagons restricted to 20km buffer
@@ -52,19 +56,24 @@ loadd(admin,
## Initial values and constants ####
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tStopSub <- 80
-initWildboarSub <- setWildBoarModelInitialValues(wbArray4sim = wildBoarArraySub,
- wildBoarObs = wildBoarObsSub,
- stepsPerChunk = stepsPerChunkSub,
- bboxAll = bboxAll) # For fitting MUST OMIT tStop ARGUMENT
-constWildboarSub <- setWildBoarModelConstants (wbArray4sim = wildBoarArraySub,
- Hex = hexSub,
- HexCentroids = hexCentroidsSub,
- Outbreaks = outbreaks,
- stepsPerChunk = stepsPerChunkSub,
- bboxAll = bboxAll,
- tStop = tStopSub,
- probNeighAS1km = probNeighAS1km,
- probNeighAS2km = probNeighAS2km)
+initWildboarSub <- setWildBoarModelInitialValues(wbArray4sim = wildBoarArraySub,
+ wildBoarObs = wildBoarObsSub,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll) # For fitting MUST OMIT tStop ARGUMENT
+constWildboarSub <- setWildBoarModelConstants (wbArray4sim = wildBoarArraySub,
+ Hex = hexSub,
+ HexCentroids = hexCentroidsSub,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll,
+ tStop = tStopSub,
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km)
+ # probNeighAS1km = probNeighAS1km,
+ # probNeighAS2km = probNeighAS2km
+
##%%%%%%%%%%%%%%%%%%%
## Build R model ####
##%%%%%%%%%%%%%%%%%%%
diff --git a/nimble/model9/model_definition.R b/nimble/model9/model_definition.R
index 31af076..42d9a2e 100644
--- a/nimble/model9/model_definition.R
+++ b/nimble/model9/model_definition.R
@@ -58,7 +58,9 @@ bugsBoarCode = nimbleCode({
pRadius[1] <- 1.0
pRadius[2] <- 0.0
radiusActiveSearch ~ dcat(pRadius[1:2]) # Can set node value to either 1 or 2 in R code, to switch between two fixed parameters
- probNeighAS <- probNeighAS1km * (radiusActiveSearch==1) + probNeighAS2km * (radiusActiveSearch==2)
+ ## probNeighAS <- probNeighAS1km * (radiusActiveSearch==1) + probNeighAS2km * (radiusActiveSearch==2)
+ weightTauAHome <- weightTauAHome1km * (radiusActiveSearch==1) + weightTauAHome2km * (radiusActiveSearch==2)
+ weightTauAAway <- weightTauAAway1km * (radiusActiveSearch==1) + weightTauAAway2km * (radiusActiveSearch==2)
## Observed &/or Simulated Data
obsY[1:(nObsStates + returnI), 1:nTimeChunksPlus1, 1:nPops] <- simulateObservations(
tIntro = tIntro, xIntro = xIntro, yIntro = yIntro,
@@ -73,14 +75,19 @@ bugsBoarCode = nimbleCode({
neighStart = neighStart[1:nPops], neighEnd = neighEnd[1:nPops],
seperatedByFence = seperatedByFence[1:lNeighVec],
omegaFence = omegaFence, tFence = tFence, lNeighVec = lNeighVec,
- inHuntZone = inHuntZone[1:nPops], probNeighAS = probNeighAS,
+ inHuntZone = inHuntZone[1:nPops],
+ # probNeighAS = probNeighAS,
+ weightTauAHome = weightTauAHome,
+ weightTauAAway = weightTauAAway,
+ radiusActiveSearch = radiusActiveSearch,
tDelayActiveSearch = tDelayActiveSearch, decayTauA = decayTauA)
# For monitor2
loglik <- 0
})
-setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks, stepsPerChunk, bboxAll, tStop, returnI=FALSE, probNeighAS1km, probNeighAS2km) {
+setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks, stepsPerChunk, bboxAll, tStop, returnI=FALSE,
+ weightTauAHome1km, weightTauAAway1km, weightTauAHome2km, weightTauAAway2km) { ## probNeighAS1km, probNeighAS2km
## wbArray4sim - an element of {wildBoarArrayAll, wildBoarArraySub}
xlims = bboxAll[c("xmin","xmax")]
ylims = bboxAll[c("ymin","ymax")]
@@ -117,8 +124,12 @@ setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks
yMax = ylims[2],
seperatedByFence = attributes(Hex)$seperatedByFence,
inHuntZone = 1.0 * Hex$inHuntZone,
- probNeighAS1km = probNeighAS1km,
- probNeighAS2km = probNeighAS2km
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km
+ ## probNeighAS1km = probNeighAS1km,
+ ## probNeighAS2km = probNeighAS2km
## x0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"X"]), # X coordinate of first observed infected farm
## y0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"Y"]) # Y coordinate of first observed infected farm
)
@@ -201,7 +212,10 @@ simulateObservations = nimbleFunction (
tFence = double(0, default=1E6),
lNeighVec = integer(0),
inHuntZone = double(1),
- probNeighAS = double(0),
+ ## probNeighAS = double(0),
+ weightTauAHome = double(0),
+ weightTauAAway = double(0),
+ radiusActiveSearch = integer(0),
tDelayActiveSearch = integer(0),
decayTauA = double(0) ) {
# Indices for states
@@ -296,7 +310,12 @@ simulateObservations = nimbleFunction (
pHome_x_pAway_x_2 = 2*pHome*pAway
pAway_x_pAway_x_2 = 2*pAway*pAway
nNeigh_x_nNeigh = nimNumeric(length=nPops, value=nNeigh[1:nPops]*nNeigh[1:nPops])
- probNeighAS_x_tauA = probNeighAS * tauA
+ ## Adjust tauA for home/away pixels and 1km/2km radius
+ ## i.e. the weighting affects the daily probability, not the rate.
+ tauAHome = -log(1 - weightTauAHome * (1 - exp(-tauA)))
+ tauAAway = -log(1 - weightTauAAway * (1 - exp(-tauA)))
+ ## probNeighAS_x_tauA = probNeighAS * tauA
+ ##
## Dynamics
for (t in floor(tIntro):tStop) {
# Indicator to turn off/on observation processes
@@ -401,10 +420,10 @@ simulateObservations = nimbleFunction (
if (indicatorUpdateAS & (dStates[Hpos] + dStates[Sa] + dStates[Sp] > 0)) {
# Then an active search will start in 3 days (delay for test) + 7 days (fixed delay for organising search in simulated data)
tauA_matrix[pop, tPlus1] = decayTauA * tauA_matrix[pop, t]
- tauA_matrix[pop, tPlusDelay] = tauA
+ tauA_matrix[pop, tPlusDelay] = tauAHome
for (ii in neighStart[pop]:neighEnd[pop]) {
iNeigh = neighVec[ii]
- tauA_matrix[iNeigh, tPlusDelay] = min(tauA, probNeighAS_x_tauA + tauA_matrix[iNeigh, tPlusDelay])
+ tauA_matrix[iNeigh, tPlusDelay] = min(tauAHome, tauAAway + tauA_matrix[iNeigh, tPlusDelay])
}
}
} # End loop on pop
diff --git a/nimble/model9/simWildBoarFromMCMC9.R b/nimble/model9/simWildBoarFromMCMC9.R
index 087628a..1c9568e 100644
--- a/nimble/model9/simWildBoarFromMCMC9.R
+++ b/nimble/model9/simWildBoarFromMCMC9.R
@@ -1,7 +1,7 @@
# source(here::here("nimble/model9/simWildBoarFromMCMC9.R"))
# rm(list=ls())
-SCENARIO <- 8
+SCENARIO <- 1
# This script jumps through MCMC output and generates weighted means & stdevs of the simulations
@@ -43,8 +43,12 @@ if (reBuildModels) {
huntZone,
hexAll,
hexCentroidsAll,
- probNeighAS1km,
- probNeighAS2km,
+ weightTauAHome1km,
+ weightTauAAway1km,
+ weightTauAHome2km,
+ weightTauAAway2km,
+ #probNeighAS1km,
+ #probNeighAS2km,
outbreaks,
pig_sites,
stepsPerChunkAll,
@@ -60,16 +64,20 @@ if (reBuildModels) {
bboxAll = bboxAll,
tStop = tStopAll,
returnI = TRUE)
- constWildboarAll = setWildBoarModelConstants(wbArray4sim = wildBoarArrayAll,
- Hex = hexAll,
- HexCentroids = hexCentroidsAll,
- Outbreaks = outbreaks,
- stepsPerChunk = stepsPerChunkAll,
- bboxAll = bboxAll,
- tStop = tStopAll,
- probNeighAS1km = probNeighAS1km,
- probNeighAS2km = probNeighAS2km,
- returnI = TRUE)
+ constWildboarAll = setWildBoarModelConstants(wbArray4sim = wildBoarArrayAll,
+ Hex = hexAll,
+ HexCentroids = hexCentroidsAll,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkAll,
+ bboxAll = bboxAll,
+ tStop = tStopAll,
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km,
+ ## probNeighAS1km = probNeighAS1km,
+ ## probNeighAS2km = probNeighAS2km,
+ returnI = TRUE)
##%%%%%%%%%%%%%%%%%%%
## Build R model ####
##%%%%%%%%%%%%%%%%%%%
@@ -216,4 +224,7 @@ for (scenario in SCENARIO) { # scenario=1
## Write scenariosList as a series of files
## load(simFile)
-## probNeighAS = probNeighAS1km
+
+
+## 02h58
+## Sys.time()
diff --git a/src/functions.R b/src/functions.R
index 208e9c1..80ede0b 100644
--- a/src/functions.R
+++ b/src/functions.R
@@ -1404,7 +1404,7 @@ makeHexTiny <- function(HexAll, CellSize) { # HexAll=hexAll; CellSize=cellSize
return(HexTiny)
}
-sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) { # HexTiny=hexTiny; CellSize=cellSize
+sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=100) { # HexTiny=hexTiny; CellSize=cellSize
radiusM <- set_units(radiusInM, m)
centroids <- st_centroid(HexTiny) %>% st_geometry()
distMat <- st_distance(centroids, centroids)
@@ -1412,9 +1412,10 @@ sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) {
HexCenter <- HexTiny[iCenter,] %>% st_geometry()
HexOuter <- st_difference(HexTiny, HexCenter) %>% select(geometry)
##
- delta <- CellSize * 0.6
+ delta <- CellSize * 0.6
##
grid = st_make_grid(HexCenter, cellsize = resolution, square = FALSE, offset = c(0, 0))
+ ## plot(HexTiny[,1])
## plot(grid, add=TRUE, col="blue")
dGrid <- data.frame(id=1:length(grid)) # data.frame(id=1:length(grid))
dGrid$geometry = grid # grid
@@ -1434,8 +1435,11 @@ sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) {
##
sumOuter <- overlapsOuter %>% st_area() %>% sum()
sumInner <- overlapsInner %>% st_area() %>% sum()
- (proportionActiveSearchInOneNeighbour <- sumOuter / sumInner)
- return(proportionActiveSearchInOneNeighbour)
+ ## Original output
+ ## (proportionActiveSearchInOneNeighbour <- sumOuter / sumInner)
+ ## return(proportionActiveSearchInOneNeighbour)
+ ## New output
+ return(c(sumInner, sumOuter))
}
diff --git a/src/plan.R b/src/plan.R
index f186c76..7074346 100644
--- a/src/plan.R
+++ b/src/plan.R
@@ -593,11 +593,18 @@ plan <- drake_plan(
## They integrate over all possible locations (within the central hexagon) for an Active Search centroid AND over all possible locations visited by the active search
## This provides a weight for the strength of an active search in a pixel that does not containing the active search centroid (i.e. a neighbour)
## These weights will be added as constants in model8 to reduce bias from the current (false) assumption that the active search is globally homogeneous
- probNeighAS1km = sampleProbActive(hexTiny, cellSize, radiusInM=1000, resolution=100)
+ innerOuterAreas1km = sampleProbActive(hexTiny, cellSize, radiusInM=1000, resolution=100)
,
- probNeighAS2km = sampleProbActive(hexTiny, cellSize, radiusInM=2000, resolution=100)
+ innerOuterAreas2km = sampleProbActive(hexTiny, cellSize, radiusInM=2000, resolution=100)
,
+ weightTauAHome1km = 1, ## innerOuterAreas1km[1]/innerOuterAreas1km[1]
+ weightTauAAway1km = innerOuterAreas1km[2] / innerOuterAreas1km[1],
+ weightTauAHome2km = innerOuterAreas2km[1] / innerOuterAreas1km[1],
+ weightTauAAway2km = innerOuterAreas2km[2] / innerOuterAreas1km[1],
+ ## probNeighAS2km =
+ ## probNeighAS1km =
+
## hexWhat = st_make_grid(hexTiny, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
## hexWhat = st_make_grid(HexCenter, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
## hexWhat <- hexWhat[1,]
--
GitLab
From 31b3d0179a1607ad2e9adf6a01f97cf78a8002c8 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 11:25:38 +0100
Subject: [PATCH 2/9] Started model10 with new active search parameterisation
in model9.
---
nimble/model10/README | 3 +
nimble/model10/fitWildBoarModel9.R | 372 +++++++++++++++++++++
nimble/model10/model_definition.R | 445 ++++++++++++++++++++++++++
nimble/model10/simWildBoarFromMCMC9.R | 230 +++++++++++++
4 files changed, 1050 insertions(+)
create mode 100644 nimble/model10/README
create mode 100644 nimble/model10/fitWildBoarModel9.R
create mode 100644 nimble/model10/model_definition.R
create mode 100644 nimble/model10/simWildBoarFromMCMC9.R
diff --git a/nimble/model10/README b/nimble/model10/README
new file mode 100644
index 0000000..44d15e1
--- /dev/null
+++ b/nimble/model10/README
@@ -0,0 +1,3 @@
+Model code adaped from model9
+
+- New parameterisation for active search.
\ No newline at end of file
diff --git a/nimble/model10/fitWildBoarModel9.R b/nimble/model10/fitWildBoarModel9.R
new file mode 100644
index 0000000..dcf581d
--- /dev/null
+++ b/nimble/model10/fitWildBoarModel9.R
@@ -0,0 +1,372 @@
+# source(here::here("nimble/model9/fitWildBoarModel9.R"))
+# rm(list=ls())
+iNSim = 4
+
+source(here::here("src/packages.R"))
+source(here::here("src/functions.R"))
+source(here::here("nimble/nimbleFunctions.R"))
+source(here::here("nimble/model9/model_definition.R"))
+options(width=1000)
+# drake::r_make() # TO UPDATE INIT AND CONST
+
+## iAttractI <- 9 # 1:9
+## (attractIVec <- ilogit(seq(-4,4,by=1)))
+## (iAttractIFixed <- attractIVec[iAttractI])
+## (logit_attractI <- logit(iAttractIFixed))
+
+## nimPrint("iAttractI=",iAttractI, " logit(attractI)=",logit_attractI)
+
+(baseDir <- here::here("nimble/model9"))
+setwd(baseDir)
+## (mcmcDir <- paste0(baseDir, "/MCMC_lpAI",logit_attractI))
+
+
+nSimVec = c(10, 50, 100, 500, 1000, 2000)
+nSim = nSimVec[iNSim]
+(mcmcDir = paste0(baseDir, paste0("/MCMC",1:length(nSimVec))[iNSim])) ## "/MCMC2"))
+
+nimPrint("nSim=",nSim)
+nimPrint("mcmcDir=",mcmcDir)
+
+if (!file.exists(mcmcDir))
+ system(paste("mkdir", mcmcDir))
+
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+## Choose which model to use - all hexagons or a subset ####
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+loadd(admin,
+ bboxAll,
+ hexSub,
+ hexAll,
+ hexCentroidsSub,
+ fence,
+ weightTauAHome1km,
+ weightTauAAway1km,
+ weightTauAHome2km,
+ weightTauAAway2km,
+ ## probNeighAS1km,
+ ## probNeighAS2km,
+ outbreaks,
+ stepsPerChunkSub,
+ wildBoarArraySub, # observations aggregated over ten day intervals, hexagons restricted to 20km buffer
+ wildBoarObsSub,
+ pig_sites)
+
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+## Initial values and constants ####
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+tStopSub <- 80
+initWildboarSub <- setWildBoarModelInitialValues(wbArray4sim = wildBoarArraySub,
+ wildBoarObs = wildBoarObsSub,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll) # For fitting MUST OMIT tStop ARGUMENT
+constWildboarSub <- setWildBoarModelConstants (wbArray4sim = wildBoarArraySub,
+ Hex = hexSub,
+ HexCentroids = hexCentroidsSub,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll,
+ tStop = tStopSub,
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km)
+ # probNeighAS1km = probNeighAS1km,
+ # probNeighAS2km = probNeighAS2km
+
+##%%%%%%%%%%%%%%%%%%%
+## Build R model ####
+##%%%%%%%%%%%%%%%%%%%
+rWildBoarModelSub = nimbleModel(bugsBoarCode,
+ constants=constWildboarSub,
+ inits=initWildboarSub,
+ calculate = FALSE, debug = FALSE) # TRUE
+
+simulate(rWildBoarModelSub, c("scTIntro", "tIntro"))
+nimPrint("Sub: tIntro = ", rWildBoarModelSub$tIntro)
+
+##%%%%%%%%%%%%%%%%%%%%%%
+## Useful node sets ####
+##%%%%%%%%%%%%%%%%%%%%%%
+detNodes = rWildBoarModelSub$getNodeNames(determOnly=TRUE)
+obsNodes = sub("\\[.*","",detNodes[grep("obsY", detNodes)])
+detNodesNoObs = detNodes[grep("obsY", detNodes, invert=TRUE)]
+stochNodes = rWildBoarModelSub$getNodeNames(stochOnly=TRUE)
+
+#################################
+## Define Target & Fixed nodes ##
+#################################
+Fixed <- c("logit_pR", "radiusActiveSearch")
+Target <- sort(c(
+ "logit_pHuntY",
+ "logit_scTIntro",
+ "logit_pIntroX",
+ "logit_pIntroY",
+ "logit_pHome",
+ "logit_attractI",
+ "logit_omegaFence",
+ "log_tauP",
+ "log_tauPArmy",
+ "log_tauA",
+ "log_tauHArmy",
+ "log_beta"
+))
+(Target)
+Monitors <- Target
+
+
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+## Initial sub-optimal initialisation ####
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (TRUE) { # FALSE
+ samplesFile <- paste0(mcmcDir,"/samples.csv")
+ if (file.exists(samplesFile)) {
+ nimPrint("Using ", samplesFile)
+ previousSamples <- read.table(file=samplesFile, header = TRUE)
+ for (iTarget in 1:length(Target)) {
+ if (is.element(Target[iTarget], colnames(previousSamples)))
+ rWildBoarModelSub[[Target[iTarget]]] <- tail(previousSamples[,Target[iTarget]],1)
+ nimPrint(Target[iTarget], " = ", rWildBoarModelSub[[Target[iTarget]]])
+ }
+ } else {
+ model1SamplesFile <- "/home/david/asf-challenge/nimble/model1/MCMC/samples.csv"
+ nimPrint("Using ", model1SamplesFile)
+ previousSamples <- read.table(file=model1SamplesFile, header = TRUE)
+ previousSamples <- tail(previousSamples,1)
+ (rWildBoarModelSub$logit_scTIntro <- logit(1-(14)/100))
+ (rWildBoarModelSub$logit_pIntroX <- previousSamples$logit_pIntroX)
+ (rWildBoarModelSub$logit_pIntroY <- previousSamples$logit_pIntroY)
+ (rWildBoarModelSub$logit_pHome <- previousSamples$logit_pHome)
+ (rWildBoarModelSub$logit_attractI <- previousSamples$logit_attractI)
+ (rWildBoarModelSub$logit_omegaFence <- logit(0.999))
+ (rWildBoarModelSub$log_tauP <- -8) # 1 + log(exp(previousSamples$log_tauCarDet) * (1-ilogit(previousSamples$logit_pActive))))
+ (rWildBoarModelSub$log_tauPArmy <- -6.6)
+ (rWildBoarModelSub$log_tauA <- 2 + log(exp(previousSamples$log_tauCarDet) * ilogit(previousSamples$logit_pActive)))
+ (rWildBoarModelSub$log_tauHArmy <- -3.47)
+ (rWildBoarModelSub$log_beta <- previousSamples$log_beta)
+ (iniPropCov <- diag(1E-3, length(Target)))
+ }
+ (rWildBoarModelSub$logit_pHuntY <- 0) # previousSamples$logit_pHuntY)
+ (rWildBoarModelSub$logit_pR <- logit(0.05))
+}
+
+
+##%%%%%%%%%%%%%%%%%%%%
+## Configure MCMC ####
+##%%%%%%%%%%%%%%%%%%%%
+(simulate(rWildBoarModelSub, detNodesNoObs))
+rWildBoarModelSub$tauH
+rWildBoarModelSub$tauHArmy
+rWildBoarModelSub$omegaFence
+rWildBoarModelSub$attractI
+
+val <- {}
+for (tt in Target) {
+ val[tt] <- eval(parse(text=paste0("rWildBoarModelSub$", tt)))
+ nimPrint(tt, " = ", val[tt])
+}
+
+##%%%%%%%%%%%%%%%%%%%%
+## Compile models ####
+##%%%%%%%%%%%%%%%%%%%%
+if (is.element("cWildBoarModelSub", ls())) {
+ nimCopy(from=rWildBoarModelSub, to=cWildBoarModelSub, nodes = stochNodes)
+} else {
+ cWildBoarModelSub = compileNimble(rWildBoarModelSub, showCompilerOutput = TRUE) # FALSE
+}
+
+#####################################################################################
+## Visual Check on Simulation - CAREFUL TO REPLACE SIMUALTED DATA WITH ORIGINAL!!! ##
+#####################################################################################
+if (TRUE) { # FALSE
+ # simulate(rWildBoarModelSub, detNodes)
+ simulate(cWildBoarModelSub, detNodes)
+ fNames <- c("1HN","2HP","3AS","4PS")
+ for (tenDayChunk in 1:9) {
+ print(data.frame(data=c(sum(wildBoarArraySub[1,tenDayChunk,]),
+ sum(wildBoarArraySub[2,tenDayChunk,]),
+ sum(wildBoarArraySub[3,tenDayChunk,]),
+ sum(wildBoarArraySub[4,tenDayChunk,])),
+ model=c(sum(cWildBoarModelSub$obsY[1,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[2,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[3,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[4,tenDayChunk,])),
+ row.names = c("HN", "HP", "AS", "PS")))
+ }
+ if (FALSE) { # TRUE
+ for (tenDayChunk in 1:8) {
+ for (iCol in 1:4) {
+ nimPrint(tenDayChunk, " ", iCol)
+ mapInfectiousDayX(t(rWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
+ Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ fileName=paste0("check_chunk", tenDayChunk, "_", fNames[iCol], "_data"),
+ StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
+ mapInfectiousDayX(t(cWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
+ Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_sim"),
+ StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
+ }
+ }
+ }
+}
+
+##%%%%%%%%%%
+## MCMC ####
+##%%%%%%%%%%
+reconfigure <- TRUE
+allowResets <- TRUE # FALSE ## Set to TRUE prior to converging, FALSE once stable
+nIter <- 1000
+for (iRun in 1:100) { # iRun=1
+ nimPrint("nSim=",nSim)
+ nimPrint("mcmcDir=",mcmcDir)
+ samplesFile <- paste0(mcmcDir,"/samples.csv")
+ samples2File <- paste0(mcmcDir,"/samples2.csv")
+ samplesExist <- file.exists(samplesFile) & file.exists(samples2File)
+ if (samplesExist) {
+ loglikSamples <- as.mcmc(read.table(file=samples2File, header = TRUE))
+ ess <- effectiveSize(loglikSamples)
+ nimPrint("Effective size of loglikSamples = ", ess)
+ if (iRun==1 | (allowResets & ess < 10)) {
+ ######################################
+ ## Verify observed data are in tact ##
+ ######################################
+ rWildBoarModelSub$obsY[,,] <- wildBoarArraySub[,,]
+ cWildBoarModelSub$obsY[,,] <- wildBoarArraySub[,,]
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Load previous samples ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ previousSamples <- read.table(file=samplesFile, header = TRUE)
+ previousSamples <- tail(previousSamples, n=round(nrow(previousSamples)/3))
+ for (iTarget in 1:length(Target)) {
+ if (is.element(Target[iTarget], colnames(previousSamples)))
+ rWildBoarModelSub[[Target[iTarget]]] <- tail(previousSamples[,Target[iTarget]],1)
+ nimPrint(Target[iTarget], " = ", rWildBoarModelSub[[Target[iTarget]]])
+ }
+ nimCopy(from=rWildBoarModelSub, to=cWildBoarModelSub, nodes = stochNodes)
+ simulate(cWildBoarModelSub, detNodesNoObs)
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Reinitialise Cov Matrix ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ iniPropCov <- cov(previousSamples)
+ diag(iniPropCov)
+ while (any(0>eigen(iniPropCov)$values)) {
+ scalePropCov <- 1.1
+ diag(iniPropCov) <- diag(iniPropCov) * scalePropCov
+ iniPropCov <- iniPropCov / scalePropCov
+ }
+ diag(iniPropCov)
+ nimPrint(iniPropCov)
+ nimPrint(cov2cor(iniPropCov))
+ reconfigure <- TRUE
+ }
+ }
+ if (iRun==1 | reconfigure) {
+ ##%%%%%%%%%%%%%%%%%%%%%%%
+ ## Reinitialise MCMC ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%
+ mcPoissonLoglik <- setupMcPoissonLoglik(model=rWildBoarModelSub, obsYNode='obsY', obsY = wildBoarArraySub,
+ constants = list(nSim=nSim, Log=TRUE), targetNodes = Target, fixedNodes = Fixed)
+ mcmcConf <- configureMCMC(rWildBoarModelSub, monitors=Monitors, monitors2="loglik", nodes = NULL)
+ mcmcConf$printSamplers()
+ mcmcConf$getMonitors()
+ mcmcConf$getMonitors2()
+ mcmcConf$removeSamplers()
+ mcmcConf$addSampler(target = Target, type = 'sampler_RW_llFunction_block_custom', ## sampler_RW_llFunction_block_custom - this gives "could not find function "extractControlElement"" which is nuts
+ control = list(llFunction = mcPoissonLoglik,
+ includesTarget = FALSE, ## FALSE tells sampler to add prior weights for target
+ propCov = iniPropCov,
+ scale = ifelse(samplesExist, 1.0, 0.01), ## Default=1 - works bad when initial covProp is terrible
+ ## adaptFactorExponent = 0.7, ## Default=0.8
+ adaptive = TRUE,
+ adaptInterval=200)) ## Default = 200
+ print(mcmcConf)
+ Cmcmc = compileNimble(buildMCMC(mcmcConf))
+ reconfigure <- FALSE
+ }
+ ##%%%%%%%%%%%%%%
+ ## Run MCMC ####
+ ##%%%%%%%%%%%%%%
+ STime <- run.time(Cmcmc$run(nIter, reset=FALSE)) ## 23h22
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Import into coda - filtering a NA line that can be generated when recompiling the modelVal##ues object ##
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ samples <- as.matrix(Cmcmc$mvSamples)
+ samples <- samples[samples[,"logit_pHome"]!=0,]
+ if(any(is.na(samples[1,])))
+ samples <- samples[-1,]
+ samples <- as.mcmc(samples)
+ ##
+ samples2 <- as.matrix(Cmcmc$mvSamples2)
+ samples2 <- samples2[samples2[,"loglik"]!=0,]
+ if(is.na(samples2[1]))
+ samples2 <- samples2[-1]
+ samples2 <- as.mcmc(samples2)
+ ##%%%%%%%%%%%%%%%%%%%%
+ ## Plot with CODA ####
+ ##%%%%%%%%%%%%%%%%%%%%
+ pdf(file=paste0(mcmcDir,"/trajectories_",nSim,".pdf"))
+ plot(samples)
+ dev.off()
+ ##
+ pdf(file=paste0(mcmcDir,"/trajectories-loglik_",nSim,".pdf"))
+ plot(samples2)
+ dev.off()
+ ##
+ write.table(samples, file=paste0(mcmcDir,"/samples.csv"), row.names = FALSE)
+ write.table(samples2, file=paste0(mcmcDir,"/samples2.csv"),row.names = FALSE)
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Check marginalised simulated values ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for (tenDayChunk in 1:9) {
+ print(data.frame(data=c(sum(wildBoarArraySub[1,tenDayChunk,]),
+ sum(wildBoarArraySub[2,tenDayChunk,]),
+ sum(wildBoarArraySub[3,tenDayChunk,]),
+ sum(wildBoarArraySub[4,tenDayChunk,])),
+ model=c(sum(cWildBoarModelSub$obsY[1,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[2,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[3,tenDayChunk,]),
+ sum(cWildBoarModelSub$obsY[4,tenDayChunk,])),
+ row.names = c("HN", "HP", "AS", "PS")))
+ }
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Check latest simulation ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ if (FALSE) {
+ for (tenDayChunk in 1:8) {
+ for (iCol in 1:4) {
+ nimPrint(tenDayChunk, " ", iCol)
+ mapInfectiousDayX(t(wildBoarArraySub[iCol,,]), Hex=hexSub,
+ Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_data"),
+ StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
+ mapInfectiousDayX(t(cWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
+ Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_sim"),
+ StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
+ }
+ }
+ }
+}
+
+
+
+
+
+### What prior can stop tauP shrinking to 0 and remain realistic ?
+## M <- 0
+## SD <- 3
+## (yMax <- qlnorm(0.99,M,SD)) ## exp(M+SD^2)
+## par(mfrow=n2mfrow(2))
+## curve(dlnorm(x,M,SD),0,yMax, n=1001)
+## (R <- exp(M - SD^2)) ## Mode of lnorm
+## ## Is mode reasonable?
+## (yMax <- qexp(0.99,R)) ## exp(M+SD^2)
+## curve(pexp(x,rate=R),0, yMax, ylim=0:1)
+## qexp(0.99,rate=R) / 365 # 102.2358 years to have 99% prob of detection
+## pgeom(yMax,pexp(1,rate=R)) # 0.9889418
+## qgeom(0.99,pexp(1,rate=R)) / 365 # 102.2356 years to have 99% prob of detection
+## log(R) # -9.21034
diff --git a/nimble/model10/model_definition.R b/nimble/model10/model_definition.R
new file mode 100644
index 0000000..42d9a2e
--- /dev/null
+++ b/nimble/model10/model_definition.R
@@ -0,0 +1,445 @@
+## TO-DO
+## - visualise relation between passive and active search
+## - adapt to hexagonal grid
+
+
+bugsCode4MV = nimbleCode({
+ logit(scTIntro) ~ dLogitBeta(7, 7) # Introduction date (on 0:1)
+ logit(pIntroX) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
+ logit(pIntroY) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
+ logit(pHome) ~ dLogitBeta(2, 2) # Connectivity: pHome=p(wild boar in home cell). (1-pHome)=p(wild boar in neighbouring cell)
+ logit(attractI) ~ dLogitBeta(2, 2) # Relative contact rate, i.e. how attractive is an infectious individual relative to a carcass.
+ logit(omegaFence) ~ dLogitBeta(2, 2) # Fence efficacy
+ log(tauP) ~ dnorm(0, sd=3) # Carcass detection rate - passive detection. TESTING THIS PRIOR TO AVOID DISAPPEARING TO ZERO
+ log(tauA) ~ dLogExp(rate = 1/1e+07) # Carcass detection rate - active detection
+ log(tauHArmy) ~ dLogExp(rate = 1/1e+07) # Increase in hunting rate around fenced area after day 60
+ log(tauPArmy) ~ dLogExp(rate = 1/1e+07) # Increase in passive search rate around fenced area after day 60
+ log(beta) ~ dLogExp(rate = 1) # The scale of transmission - multiply with the realtive attractivity of Infectious and Carcasses.
+ logit(pHuntY) ~ dLogitBeta(330, 330)
+})
+
+bugsBoarCode = nimbleCode({
+ ## Known Parameters
+ omega <- 0.8 # Habitat preference
+ tauI <- 1/7 # Incubation rate
+ tauC <- 1/7 # Disease induced mortality rate
+ tauR <- tauC * pR/(1 - pR) # Recovery rate so that pR (recovfery probability ) is 0.05
+ tauRot <- 1/90 # Carcass decomposition rate
+ tauH <- -log((1 - pHuntY)^(12/(8 * 365))) # Hunting rate. Hunters think pHuntY=0.5 (prob. shot in one year). Hunting season is 8 month long.
+ pTest <- 20/100 # Probability to test a hunted wild boar
+ tDelayActiveSearch <- 10 # 3 days for the test and 7 days to organise search. 100% constant in their simulaed data.
+ decayTauA <- (1 - 1/7) # Geometric decay of tauA
+ ## Unknown Parameters
+ logit(scTIntro) ~ dLogitBeta(7, 7) # Introduction date (on 0:1)
+ logit(pIntroX) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
+ logit(pIntroY) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
+# logit(pHome) ~ dLogitBeta(2, 2) # Connectivity: pHome=p(wild boar in home cell). (1-pHome)=p(wild boar in neighbouring cell)
+ logit(attractI) ~ dLogitBeta(2, 2) # Relative contact rate, i.e. how attractive is an infectious individual relative to a carcass.
+ logit(omegaFence) ~ dLogitBeta(2, 2) # Fence efficacy
+ log(tauP) ~ dnorm(0, sd=3) # Carcass detection rate - passive detection. TESTING THIS PRIOR TO AVOID DISAPPEARING TO ZERO
+ log(tauA) ~ dLogExp(rate = 1/1e+07) # Carcass detection rate - active detection
+ log(tauHArmy) ~ dLogExp(rate = 1/1e+07) # Increase in hunting rate around fenced area after day 60
+ log(tauPArmy) ~ dLogExp(rate = 1/1e+07) # Increase in passive search rate around fenced area after day 60
+# log(beta) ~ dLogExp(rate = 1) # The scale of transmission - multiply with the realtive attractivity of Infectious and Carcasses.
+ # Changes applied to MCMC4 15h Sunday to attempt reducing beta~pHome redundancy
+ log(beta) ~ dLogExp(rate = 0.1)
+ logit(pHome) ~ dLogitBeta(10, 10)
+ ##
+ betaI <- beta * attractI
+ betaC <- beta * (1 - attractI)
+ tIntro <- tMin - scTIntro * tMin
+ xIntro <- xMin + pIntroX * (xMax - xMin)
+ yIntro <- yMin + pIntroY * (yMax - yMin)
+ ## Somewhat known, can be fixed for simplicity
+ logit(pR) ~ dLogitBeta(5/2.18, 95/2.18)
+ logit(pHuntY) ~ dLogitBeta(330, 330)
+ ## Other constants
+ tFence <- 60 # Last day without fence and extra hunting
+ pRadius[1] <- 1.0
+ pRadius[2] <- 0.0
+ radiusActiveSearch ~ dcat(pRadius[1:2]) # Can set node value to either 1 or 2 in R code, to switch between two fixed parameters
+ ## probNeighAS <- probNeighAS1km * (radiusActiveSearch==1) + probNeighAS2km * (radiusActiveSearch==2)
+ weightTauAHome <- weightTauAHome1km * (radiusActiveSearch==1) + weightTauAHome2km * (radiusActiveSearch==2)
+ weightTauAAway <- weightTauAAway1km * (radiusActiveSearch==1) + weightTauAAway2km * (radiusActiveSearch==2)
+ ## Observed &/or Simulated Data
+ obsY[1:(nObsStates + returnI), 1:nTimeChunksPlus1, 1:nPops] <- simulateObservations(
+ tIntro = tIntro, xIntro = xIntro, yIntro = yIntro,
+ nObsStates = nObsStates, returnI = returnI, nPops = nPops, nTimeChunksPlus1 = nTimeChunksPlus1,
+ stepsPerChunk = stepsPerChunk, xCentroid = xCentroid[1:nPops],
+ yCentroid = yCentroid[1:nPops], weightF = weightF[1:nPops],
+ weightA = weightA[1:nPops], areaHuntBag = areaHuntBag[1:nPops],
+ omega = omega, betaI = betaI, betaC = betaC, pHome = pHome,
+ pR = pR, tauA = tauA, tauP = tauP, tauPArmy = tauPArmy, tauI = tauI, tauC = tauC,
+ tauR = tauR, tauRot = tauRot, tauH = tauH, tauHArmy = tauHArmy, pHuntY = pHuntY,
+ pTest = pTest, nNeigh = nNeigh[1:nPops], neighVec = neighVec[1:lNeighVec],
+ neighStart = neighStart[1:nPops], neighEnd = neighEnd[1:nPops],
+ seperatedByFence = seperatedByFence[1:lNeighVec],
+ omegaFence = omegaFence, tFence = tFence, lNeighVec = lNeighVec,
+ inHuntZone = inHuntZone[1:nPops],
+ # probNeighAS = probNeighAS,
+ weightTauAHome = weightTauAHome,
+ weightTauAAway = weightTauAAway,
+ radiusActiveSearch = radiusActiveSearch,
+ tDelayActiveSearch = tDelayActiveSearch, decayTauA = decayTauA)
+ # For monitor2
+ loglik <- 0
+})
+
+
+setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks, stepsPerChunk, bboxAll, tStop, returnI=FALSE,
+ weightTauAHome1km, weightTauAAway1km, weightTauAHome2km, weightTauAAway2km) { ## probNeighAS1km, probNeighAS2km
+ ## wbArray4sim - an element of {wildBoarArrayAll, wildBoarArraySub}
+ xlims = bboxAll[c("xmin","xmax")]
+ ylims = bboxAll[c("ymin","ymax")]
+ nObsStates = dim(wbArray4sim)[1]
+ if (missing(tStop)) {
+ nTimeChunksPlus1 = dim(wbArray4sim)[2]
+ } else {
+ nTimeChunks = ceiling(tStop / stepsPerChunk)
+ nTimeChunksPlus1 = nTimeChunks + 1
+ }
+ nPops = dim(wbArray4sim)[3]
+ if (nPops != nrow(Hex)) stop("The following must be true: dim(wbArray4sim)[3] == nrow(Hex)")
+ #
+ const = list(
+ nPops = nPops, # Number of hexagonal patches
+ nObsStates = nObsStates, # Number of observation types for wild boar, "NT" "PT" "AS" "PS" (we drop NAS)
+ returnI = returnI, # Indicator: FALSE - simulations don't return I, TRUE - simulations do return I
+ nTimeChunksPlus1 = nTimeChunksPlus1, # Number of bins for aggregating time-step data
+ stepsPerChunk = stepsPerChunk, # Number of timesteps for each aggregated bin
+ weightF = Hex$weightF, # Weight for estimating wild boar density (from admin area's hunting bag) based on forest cover in hexagon / total forest cover in admin area.
+ weightA = Hex$weightA, # Weight for estimating wild boar density (from admin area's hunting bag) based on agro cover in hexagon / total agro cover in admin area.
+ areaHuntBag = Hex$adminAreaHuntBag, # Hunting bag data - to be divided by all hexagons from a given admin area
+ nNeigh = Hex$nNeigh, # Number of neighbours for each hexagon
+ neighStart = Hex$neighStart, # Indicates start of sub-vector of neighVec for each hexagonal patch
+ neighEnd = Hex$neighEnd, # Indicates end of sub-vector of neighVec for each hexagonal patch
+ neighVec = attributes(Hex)$neighVec, # A vector containing indices for neighbours
+ lNeighVec = length(attributes(Hex)$neighVec),
+ xCentroid = (HexCentroids %>% st_coordinates())[,"X"], # X component of each hexagon's centroid
+ yCentroid = (HexCentroids %>% st_coordinates())[,"Y"], # Y component of each hexagon's centroid
+ tMin = -100, # Lower bound on carcass-fall-from-space day
+ xMin = xlims[1],
+ xMax = xlims[2],
+ yMin = ylims[1],
+ yMax = ylims[2],
+ seperatedByFence = attributes(Hex)$seperatedByFence,
+ inHuntZone = 1.0 * Hex$inHuntZone,
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km
+ ## probNeighAS1km = probNeighAS1km,
+ ## probNeighAS2km = probNeighAS2km
+ ## x0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"X"]), # X coordinate of first observed infected farm
+ ## y0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"Y"]) # Y coordinate of first observed infected farm
+ )
+ const
+}
+
+setWildBoarModelInitialValues <- function (wbArray4sim, wildBoarObs, stepsPerChunk, tStop, bboxAll, returnI=FALSE) {
+ ## browser()
+ ## wbArray4sim - an element of {wildBoarArrayAll, wildBoarArraySub}
+ nObsStates = dim(wbArray4sim)[1]
+ nTimeChunksPlus1 = dim(wbArray4sim)[2]
+ nPops = dim(wbArray4sim)[3] # = nrow(Hex) # Same value
+ nTimeChunks = nTimeChunksPlus1 - 1
+ nReturnStates = nObsStates + returnI
+ #
+ if (missing(tStop)) {
+ nimPrint("USING ORIGINAL OBSERVATIONS IN setWildBoarModelInitialValues")
+ nimPrint("dim(obsY)[2] = ", nTimeChunksPlus1)
+ tStop = nTimeChunks * stepsPerChunk
+ myObsY = nimArray(value=0, dim=c(nReturnStates, nTimeChunksPlus1, nPops))
+ myObsY[1:nObsStates, 1:nTimeChunksPlus1, 1:nPops] = wbArray4sim[1:nObsStates, 1:nTimeChunksPlus1, 1:nPops] # These nodes are first initialised with observed data, but the data will be replaced each time this node is simulated.
+ } else if (nTimeChunks <= ceiling(tStop / stepsPerChunk)) {
+ nTimeChunks = ceiling(tStop / stepsPerChunk)
+ nTimeChunksPlus1 = nTimeChunks + 1
+ nimPrint("NOT USING ORIGINAL OBSERVATIONS IN setWildBoarModelInitialValues. IF THIS IS NOT WHAT YOU WANT THEN DELETE tStop ARGUMENT.")
+ nimPrint("setting dim(obsY)[2] to ", nTimeChunksPlus1)
+ myObsY = nimArray(value=0, dim=c(nReturnStates, nTimeChunksPlus1, nPops))
+ } else {
+ stop("Please ensure ceiling(tStop / stepsPerChunk) >= dim(wbArray4sim)[2] in order to simulate at least as far as original data.")
+ }
+ # Ball-park Initial value for location of first infectious carcass
+ xIntro = (wildBoarObs %>% filter(is.element(DET, c("AS","PS","PT"))) %>% st_coordinates() %>% colMeans())["X"]
+ yIntro = (wildBoarObs %>% filter(is.element(DET, c("AS","PS","PT"))) %>% st_coordinates() %>% colMeans())["Y"]
+ xMin = bboxAll["xmin"]
+ xMax = bboxAll["xmax"]
+ yMin = bboxAll["ymin"]
+ yMax = bboxAll["ymax"]
+ #
+ inits = list(
+ obsY = myObsY[1:nReturnStates, 1:nTimeChunksPlus1, 1:nPops],
+ logit_scTIntro = logit(0.75), # Time first infectious carcass fell from space
+ logit_pIntroX = logit((xIntro-xMin)/(xMax-xMin)),
+ logit_pIntroY = logit((yIntro-yMin)/(yMax-yMin)),
+ logit_pHome = logit(1/7), # Connectivity
+ logit_attractI = logit(0.5), # Attraction of Infectious indivuidual (relative to a Carcass)
+ logit_omegaFence = logit(0.999),
+ logit_pHuntY = logit(0.5), # Probability a wild boar is shot within one year's 8 month hunting season
+ logit_pR = logit(0.05), # Recovery probability. FIXED
+ log_beta = -3.5, # Transmission coeficient
+ log_tauP = log(0.0001), # Detection rate for passive detection
+ log_tauA = log(0.001), # Detection rate for active detection
+ log_tauHArmy = -7,
+ log_tauPArmy = -7,
+ radiusActiveSearch = 1
+ )
+ inits
+}
+
+
+######################################################
+## Stochastic simulation functions - wild boar only ##
+######################################################
+simulateObservations = nimbleFunction (
+ run = function (tIntro = double(0), xIntro = double(0), yIntro = double(0),
+ nObsStates = integer(0),
+ returnI = logical(0), ## If TRUE, an additional output for I is returned.
+ nPops = integer(0), nTimeChunksPlus1 = double(0), stepsPerChunk = double(0),
+ xCentroid = double(1), yCentroid = double(1),
+ weightF = double(1), weightA = double(1),
+ areaHuntBag = double(1), omega = double(0), betaI = double(0), betaC = double(0),
+ pHome = double(0), pR = double(0), tauA = double(0),
+ tauP = double(0), tauPArmy = double(0),
+ tauI = double(0), tauC = double(0), tauR = double(0),
+ tauRot = double(0), tauH = double(0),
+ tauHArmy = double(0), pHuntY = double(0), pTest = double(0),
+ nNeigh = double(1), neighVec = double(1),
+ neighStart = double(1), neighEnd = double(1),
+ seperatedByFence = double(1), # Must be same length as neighVec
+ omegaFence = double(0),
+ tFence = double(0, default=1E6),
+ lNeighVec = integer(0),
+ inHuntZone = double(1),
+ ## probNeighAS = double(0),
+ weightTauAHome = double(0),
+ weightTauAAway = double(0),
+ radiusActiveSearch = integer(0),
+ tDelayActiveSearch = integer(0),
+ decayTauA = double(0) ) {
+ # Indices for states
+ S = 1
+ E = 2
+ I = 3
+ R = 4
+ C = 5
+ Rot = 6
+ Hu = 7
+ Hneg = 8 # NT in docs
+ Hpos = 9 # PT in docs
+ Sa = 10 # AS in docs
+ Sp = 11 # PS in docs
+ nStates = 11
+ obState1 = nStates - nObsStates + 1
+ obState4 = nStates
+ zero4Obs = nimNumeric(value=0, length=nObsStates)
+ nReturnStates = nObsStates + returnI
+ # Indices for rates - individual
+ S_Infection = 1
+ S_HuntTestNeg = 2
+ S_HuntUntest = 3
+ E_Incubation = 4
+ E_HuntTestPos = 5
+ E_HuntUntest = 6
+ I_Recover = 7
+ I_HuntTestPos = 8
+ I_HuntUntest = 9
+ I_Death = 10
+ C_Rot = 11
+ C_DetActive = 12
+ C_DetPassive = 13
+ nRates = 13
+ # Indices for rates - vector
+ iSrates = 1:3
+ iErates = 4:6
+ iIrates = 7:10
+ iCrates = 11:13
+ iSmin = 1
+ iSmax = 3
+ iEmin = 4
+ iEmax = 6
+ iImin = 7
+ iImax = 10
+ iCmin = 11
+ iCmax = 13
+ ## Establish
+ dtDynamics1 = 1.0 # Resolution of discretisation of dynamic process
+ dtDynamics = abs(ceiling(tIntro) - tIntro) # Duration of process in first time step only
+ if (dtDynamics==0)
+ dtDynamics = dtDynamics1
+ nTimeChunks = nTimeChunksPlus1 - 1
+ tStop = nTimeChunks * stepsPerChunk * dtDynamics1 # stepsPerChunk is nb. time-steps per chunk
+ tFence1 = round(tFence + 1) ## Fence takes effect on day after completion date.
+ tStopMinusDelay = tStop - tDelayActiveSearch
+ ## breaks = nimInteger(length=nTimeChunks)
+ ## for (br in 1:nTimeChunks)
+ ## breaks[br] = br * stepsPerChunk
+ ## breaks = c(0, breaks)
+ breaks = nimInteger(length=nTimeChunksPlus1)
+ for (br in 1:nTimeChunksPlus1)
+ breaks[br] = (br-1) * stepsPerChunk
+ ## Identify location (hexagon) of initial carcass
+ distIntro = nimNumeric(length=nPops)
+ for (pop in 1:nPops) {
+ distIntro[pop] = sqrt((xIntro-xCentroid[pop])^2 + (yIntro-yCentroid[pop])^2)
+ }
+ hexIntro = which(distIntro[1:nPops]==min(distIntro[1:nPops])) ### nimPrint("hexIntro=",hexIntro)
+ ## Initialise arrays for dynamics
+ obs = nimArray(value=0, dim=c(nReturnStates, nTimeChunksPlus1, nPops))
+ states = nimMatrix(value=0, nrow=nStates, ncol=nPops)
+ pStates = nimMatrix(value=0, nrow=nStates, ncol=nPops)
+ tauA_matrix = nimMatrix(value=0, nrow=nPops, ncol=tStop) ### Index columns with max(t, 1) to acount for t<1
+ dStates = nimNumeric(value=0, length=nStates)
+ rates = nimNumeric(value=0, length=nRates)
+ vecLeave = nimNumeric(value=0, length=nRates)
+ tauA_local = 0
+ states[C, hexIntro] = 1
+ fenceInducedConnectivityReduction = nimNumeric(value=0, length=lNeighVec)
+ ## Estimate expected boar density & sample
+ ## browser()
+ Eboar = nimNumeric(value=0, length=nPops)
+ Eboar[1:nPops] = (omega*weightF[1:nPops] + (1-omega)*weightA[1:nPops]) * areaHuntBag[1:nPops] ## / pHuntY # E[boar] * pHuntY = annual_hunting_bag
+ for (pop in 1:nPops) {
+ states[S, pop] = rpois(n=1, lambda=Eboar[pop])
+ }
+ ## Weights for being at home or at a neighbour's
+ pAway = 1-pHome
+ pHome_x_pHome = pHome*pHome
+ pAway_x_pAway = pAway*pAway
+ pHome_x_pAway_x_2 = 2*pHome*pAway
+ pAway_x_pAway_x_2 = 2*pAway*pAway
+ nNeigh_x_nNeigh = nimNumeric(length=nPops, value=nNeigh[1:nPops]*nNeigh[1:nPops])
+ ## Adjust tauA for home/away pixels and 1km/2km radius
+ ## i.e. the weighting affects the daily probability, not the rate.
+ tauAHome = -log(1 - weightTauAHome * (1 - exp(-tauA)))
+ tauAAway = -log(1 - weightTauAAway * (1 - exp(-tauA)))
+ ## probNeighAS_x_tauA = probNeighAS * tauA
+ ##
+ ## Dynamics
+ for (t in floor(tIntro):tStop) {
+ # Indicator to turn off/on observation processes
+ indicatorAfterDay0 = (t > 0)
+ indicatorDay0OrAfter = (t >= 0)
+ indicatorAfterDay60 = (t > tFence)
+ indicatorUpdateAS = indicatorAfterDay0 & (t <= tStopMinusDelay) # Can do an active search before tStop+1
+ tPlusDelay = t + tDelayActiveSearch
+ tPlus1 = t + 1
+ # Vector to scale (reduce) connectivity at fence
+ if (t == tFence1)
+ fenceInducedConnectivityReduction[1:lNeighVec] = (1 - seperatedByFence[1:lNeighVec] * omegaFence)
+ # Save state of system as in previous time step
+ pStates[1:nStates, 1:nPops] = states[1:nStates, 1:nPops] # p = previous
+ # nimPrint(sum(pStates[S, 1:nPops]), " ", sum(pStates[Hneg, 1:nPops]))
+ #
+ for (pop in 1:nPops) {
+ ## Share-cell (i.e. "contact") probabilities (determined by epsillon) between i (pop) and j (neighbour)
+ pContactIi = pHome_x_pHome + pAway_x_pAway/nNeigh[pop]
+ pContactIj = pHome_x_pAway_x_2/nNeigh[pop] + 2*pAway_x_pAway_x_2 / (nNeigh_x_nNeigh[pop]) ## Assumes all pop's neighbours have same number of neighbours as pop. This is fine if hexagons overlap the sea, and thus have zero wild boar.
+ pContactCi = pHome
+ pContactCj = pAway/nNeigh[pop]
+ ## Number of contagious neighbours
+ if ( !indicatorAfterDay60 ) { ## (t <= tFence) {
+ neighSumI = sum(pStates[I, neighVec[neighStart[pop]:neighEnd[pop]]])
+ neighSumC = sum(pStates[C, neighVec[neighStart[pop]:neighEnd[pop]]])
+ } else {
+ neighSumI = sum(pStates[I, neighVec[neighStart[pop]:neighEnd[pop]]] * fenceInducedConnectivityReduction[neighStart[pop]:neighEnd[pop]])
+ neighSumC = sum(pStates[C, neighVec[neighStart[pop]:neighEnd[pop]]] * fenceInducedConnectivityReduction[neighStart[pop]:neighEnd[pop]])
+ }
+ ## Set "local" versions of rates. These should account for space-time differences in things like hunting pressure, active search prob etc etc
+ flagExtraHuntingAndTesting = indicatorAfterDay60 * inHuntZone[pop]
+ additionalHuntingPressure = flagExtraHuntingAndTesting * tauHArmy ## Will be zero unless in extra hunt zone and after day 60
+ additionalPassiveSearch = flagExtraHuntingAndTesting * tauPArmy ## Will be zero unless in extra hunt zone and after day 60
+ tauH_local = tauH + additionalHuntingPressure
+ pTest_local = pTest * indicatorAfterDay0
+ if (flagExtraHuntingAndTesting == 1)
+ pTest_local = 1.0
+ if (indicatorAfterDay0)
+ tauA_local = tauA_matrix[pop, t] * (1-flagExtraHuntingAndTesting) # tauA
+ tauP_local = (tauP + additionalPassiveSearch) * indicatorAfterDay0
+ ## Reaction Rates
+ rates[S_Infection] = pStates[S,pop] * ( betaI*(pContactIi*pStates[I,pop] + pContactIj*neighSumI) + betaC*(pContactCi*pStates[C,pop] + pContactCj*neighSumC) )
+ rates[E_Incubation] = pStates[E,pop] * tauI
+ rates[I_Recover] = pStates[I,pop] * tauR
+ rates[I_Death] = pStates[I,pop] * tauC
+ rates[C_Rot] = indicatorDay0OrAfter * pStates[C,pop] * tauRot ## The flag should prevent the initial carcass vanishing before infecting any other animal
+ rates[S_HuntUntest] = pStates[S,pop] * tauH_local * (1-pTest_local)
+ rates[E_HuntUntest] = pStates[E,pop] * tauH_local * (1-pTest_local)
+ rates[I_HuntUntest] = pStates[I,pop] * tauH_local * (1-pTest_local)
+ rates[S_HuntTestNeg] = pStates[S,pop] * tauH_local * pTest_local
+ rates[E_HuntTestPos] = pStates[E,pop] * tauH_local * pTest_local
+ rates[I_HuntTestPos] = pStates[I,pop] * tauH_local * pTest_local
+ rates[C_DetActive] = pStates[C,pop] * tauA_local
+ rates[C_DetPassive] = pStates[C,pop] * tauP_local
+ # Reactions - Part 1
+ nLeaveS = rpois(n=1, lambda = dtDynamics * (rates[S_Infection] + rates[S_HuntTestNeg] + rates[S_HuntUntest]))
+ nLeaveE = rpois(n=1, lambda = dtDynamics * (rates[E_Incubation] + rates[E_HuntTestPos] + rates[E_HuntUntest]))
+ nLeaveI = rpois(n=1, lambda = dtDynamics * (rates[I_Recover] + rates[I_HuntTestPos] + rates[I_HuntUntest]) + rates[I_Death])
+ nLeaveC = rpois(n=1, lambda = dtDynamics * (rates[C_Rot] + rates[C_DetActive] + rates[C_DetPassive]))
+ nLeaveS = min(nLeaveS, pStates[S,pop])
+ nLeaveE = min(nLeaveE, pStates[E,pop])
+ nLeaveI = min(nLeaveI, pStates[I,pop])
+ nLeaveC = min(nLeaveC, pStates[C,pop])
+ ## Reactions - Part 2
+ if (nLeaveS > 0) {
+ vecLeave[iSmin:iSmax] = rmulti(1, size=nLeaveS, prob=rates[iSmin:iSmax])
+ } else {
+ vecLeave[iSmin:iSmax] = 0 * vecLeave[iSmin:iSmax]
+ }
+ if (nLeaveE > 0) {
+ vecLeave[iEmin:iEmax] = rmulti(n=1, size=nLeaveE, prob=rates[iEmin:iEmax])
+ } else {
+ vecLeave[iEmin:iEmax] = (0 * vecLeave[iEmin:iEmax])
+ }
+ if (nLeaveI > 0) {
+ vecLeave[iImin:iImax] = rmulti(n=1, size=nLeaveI, prob=rates[iImin:iImax])
+ } else {
+ vecLeave[iImin:iImax] = (0 * vecLeave[iImin:iImax])
+ }
+ if (nLeaveC > 0) {
+ vecLeave[iCmin:iCmax] = rmulti(n=1, size=nLeaveC, prob=rates[iCmin:iCmax])
+ } else {
+ vecLeave[iCmin:iCmax] = (0 * vecLeave[iCmin:iCmax])
+ }
+ ## Reactions - Part 3
+ dStates[S] = ( - sum(vecLeave[iSmin:iSmax]) )
+ dStates[E] = (vecLeave[S_Infection] - sum(vecLeave[iEmin:iEmax]) )
+ dStates[I] = (vecLeave[E_Incubation] - sum(vecLeave[iImin:iImax]) )
+ dStates[R] = (vecLeave[I_Recover] )
+ dStates[C] = (vecLeave[I_Death] - sum(vecLeave[iCmin:iCmax]))
+ dStates[Rot] = (vecLeave[C_Rot] )
+ dStates[Hneg] = (vecLeave[S_HuntTestNeg] )
+ dStates[Hpos] = (vecLeave[I_HuntTestPos] + vecLeave[E_HuntTestPos])
+ dStates[Hu] = (vecLeave[S_HuntUntest] + vecLeave[E_HuntUntest] + vecLeave[I_HuntUntest])
+ dStates[Sa] = (vecLeave[C_DetActive] )
+ dStates[Sp] = (vecLeave[C_DetPassive] )
+ ## Update states
+ pStates[(nStates-nObsStates)+(1:nObsStates),pop] = zero4Obs[1:nObsStates] ## Set observation states to zero so that they do not accumulate in next step
+ states[1:nStates,pop] = pStates[1:nStates,pop] + dStates[1:nStates]
+ ## Update tauA_matrix
+ if (indicatorUpdateAS & (dStates[Hpos] + dStates[Sa] + dStates[Sp] > 0)) {
+ # Then an active search will start in 3 days (delay for test) + 7 days (fixed delay for organising search in simulated data)
+ tauA_matrix[pop, tPlus1] = decayTauA * tauA_matrix[pop, t]
+ tauA_matrix[pop, tPlusDelay] = tauAHome
+ for (ii in neighStart[pop]:neighEnd[pop]) {
+ iNeigh = neighVec[ii]
+ tauA_matrix[iNeigh, tPlusDelay] = min(tauAHome, tauAAway + tauA_matrix[iNeigh, tPlusDelay])
+ }
+ }
+ } # End loop on pop
+ if (indicatorDay0OrAfter) {
+ iChunk = min(which(t<=breaks))
+ obs[1:nObsStates, iChunk, 1:nPops] = obs[1:nObsStates, iChunk, 1:nPops] + states[(nStates-nObsStates)+(1:nObsStates), 1:nPops]
+ if (returnI) {
+ obs[nReturnStates, iChunk, 1:nPops] = obs[nReturnStates, iChunk, 1:nPops] + states[I, 1:nPops]
+ }
+ }
+ dtDynamics = dtDynamics1 # Sets dt for all subsequent time steps
+ } # End loop on t
+ if (returnI) {
+ obs[nReturnStates, 1:nTimeChunksPlus1, 1:nPops] = obs[nReturnStates, 1:nTimeChunksPlus1, 1:nPops] / stepsPerChunk ## Return expected value of I over each chunk
+ }
+ returnType(double(3))
+ return(obs[1:nReturnStates, 1:nTimeChunksPlus1, 1:nPops])
+ }
+)
diff --git a/nimble/model10/simWildBoarFromMCMC9.R b/nimble/model10/simWildBoarFromMCMC9.R
new file mode 100644
index 0000000..1c9568e
--- /dev/null
+++ b/nimble/model10/simWildBoarFromMCMC9.R
@@ -0,0 +1,230 @@
+# source(here::here("nimble/model9/simWildBoarFromMCMC9.R"))
+# rm(list=ls())
+
+SCENARIO <- 1
+
+# This script jumps through MCMC output and generates weighted means & stdevs of the simulations
+
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+## Set directories and MCMC output file ####
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+library(here)
+## help.start()
+(modelDir = here::here("nimble/model9")) # For selecting model version
+(mcmcDir = paste0(modelDir,"/MCMC4")) # 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,"/sims")) # For selecting model version
+setwd(modelDir)
+# setwd(mcmcDir)
+mcmcFile = paste0(mcmcDir,"/samples.csv") ## Assumed to be parameters
+mcmcFile2 = paste0(mcmcDir,"/samples2.csv") ## Assumed to be posterior log-likelihoods
+## modelInitialisationFile = paste0("../", dir("../")[grep("initialisation_model", dir("../"))])
+
+if (!file.exists(simDir))
+ system(paste("mkdir", simDir))
+
+radiusActiveSearch = 1 # 1
+
+##%%%%%%%%%%%%%%%%%%%%%%
+## Build the models ####
+##%%%%%%%%%%%%%%%%%%%%%%
+reBuildModels = TRUE # FALSE
+if (reBuildModels) {
+ setwd(modelDir)
+ source(here::here("src/packages.R"))
+ source(here::here("src/functions.R"))
+ source(here::here("nimble/nimbleFunctions.R"))
+ source(paste0(modelDir, "/model_definition.R")) # THIS ONLY LOADS FUNCTIONS (WRITEN FOR DRAKE), NOT THE ACTUAL VALUES
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Choose which model to use - all hexagons or a subset ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ loadd(admin,
+ bboxAll,
+ fence,
+ huntZone,
+ hexAll,
+ hexCentroidsAll,
+ weightTauAHome1km,
+ weightTauAAway1km,
+ weightTauAHome2km,
+ weightTauAAway2km,
+ #probNeighAS1km,
+ #probNeighAS2km,
+ outbreaks,
+ pig_sites,
+ stepsPerChunkAll,
+ wildBoarArrayAll, # observations aggregated over one day intervals, hexagons cover whole island
+ wildBoarObsAll)
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Initial values and constants ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ tStopAll = 140 ## For 30 day prediction beyond day 80
+ initWildboarAll = setWildBoarModelInitialValues(wbArray4sim = wildBoarArrayAll,
+ wildBoarObs = wildBoarObsAll,
+ stepsPerChunk = stepsPerChunkAll,
+ bboxAll = bboxAll,
+ tStop = tStopAll,
+ returnI = TRUE)
+ constWildboarAll = setWildBoarModelConstants(wbArray4sim = wildBoarArrayAll,
+ Hex = hexAll,
+ HexCentroids = hexCentroidsAll,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkAll,
+ bboxAll = bboxAll,
+ tStop = tStopAll,
+ weightTauAHome1km = weightTauAHome1km,
+ weightTauAAway1km = weightTauAAway1km,
+ weightTauAHome2km = weightTauAHome2km,
+ weightTauAAway2km = weightTauAAway2km,
+ ## probNeighAS1km = probNeighAS1km,
+ ## probNeighAS2km = probNeighAS2km,
+ returnI = TRUE)
+ ##%%%%%%%%%%%%%%%%%%%
+ ## Build R model ####
+ ##%%%%%%%%%%%%%%%%%%%
+ rWildBoarModelAll = nimbleModel(bugsBoarCode,
+ constants=constWildboarAll,
+ inits=initWildboarAll,
+ calculate = FALSE, debug = FALSE)
+ rModel4MV = nimbleModel(bugsCode4MV) ## So we can make a modelValues without the huge obsY array
+ ##%%%%%%%%%%%%%%%%%%%%%%
+ ## Useful node sets ####
+ ##%%%%%%%%%%%%%%%%%%%%%%
+ detNodes = rWildBoarModelAll$getNodeNames(determOnly=TRUE)
+ obsNodes = sub("\\[.*","",detNodes[grep("obsY", detNodes)])
+ detNodesNoObs = detNodes[grep("obsY", detNodes, invert=TRUE)]
+ stochNodes = rWildBoarModelAll$getNodeNames(stochOnly=TRUE)
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Ensure determined Nodes are initialised ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ simulate(rWildBoarModelAll, detNodesNoObs)
+ for (dn in 1:length(detNodesNoObs))
+ nimPrint(detNodesNoObs[dn], " ", rWildBoarModelAll[[detNodesNoObs[dn]]])
+ for (sn in 1:length(stochNodes))
+ nimPrint(stochNodes[sn], " ", rWildBoarModelAll[[stochNodes[sn]]])
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Compile model & initialise deterministic nodes ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ cWildBoarModelAll = compileNimble(rWildBoarModelAll, showCompilerOutput = TRUE) # FALSE
+ cWildBoarModelAll$radiusActiveSearch <- radiusActiveSearch
+ simulate(cWildBoarModelAll, detNodesNoObs)
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## Set Dimension Constants ####
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ nReturnTypes = dim(cWildBoarModelAll$obsY)[1]
+ nObsDays = dim(wildBoarArrayAll)[2]
+ nPops = dim(wildBoarArrayAll)[3]
+ nSimDays = dim(cWildBoarModelAll$obsY)[2]
+ nObsStates = nReturnTypes - constWildboarAll$returnI
+ iInfectious = nReturnTypes # Index for infectious wild boar in the output array
+}
+
+
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+## Load MCMC mcmcSamples and (optionally) remove burnin ####
+##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+mcmcSamples = read.table(file=mcmcFile, header = TRUE)
+mcmcSamples2 = read.table(file=mcmcFile2, header = TRUE) # Assumed to be posterior likelihood
+if (FALSE) { # TRUE
+ (burnin <- firstSampleAboveMean(as.numeric(mcmcSamples2[,1])))
+ mcmcSamples <- mcmcSamples[-(1:burnin),]
+ mcmcSamples2 <- mcmcSamples2[-(1:burnin),]
+}
+nimPrint("Effective sample size: ", effectiveSize(mcmcSamples))
+if (FALSE) { # TRUE
+ trajectoriesPDF = TRUE # FALSE
+ # mcmcSamples = mcmcSamples[!mcmcSamples$log_beta==0,] ## Filters out non-samples where I stopped the MCMC
+ if (trajectoriesPDF)
+ pdf(file=paste0(mcmcDir, "/trajectories.pdf"))
+ plot(as.mcmc(mcmcSamples))
+ if (trajectoriesPDF)
+ dev.off()
+ dim(mcmcSamples)
+}
+
+##%%%%%%%%%%%%%%%%%%%%%
+## Setup simulator ####
+##%%%%%%%%%%%%%%%%%%%%%
+rSimFromMCMC <- setupSimFromMCMC(model = rWildBoarModelAll,
+ model4MV = rModel4MV,
+ obsYNode = "obsY",
+ obsY = wildBoarArrayAll,
+ mcmcSamples = mcmcSamples)
+cSimFromMCMC <- compileNimble(rSimFromMCMC)
+
+## To recover original data in obsY
+## nimCopy(from=rWildBoarModelAll, to=cWildBoarModelAll, nodes="obsY")
+
+##%%%%%%%%%%%%%%%%%%%%%
+## Run simulations ####
+##%%%%%%%%%%%%%%%%%%%%%
+nSim = 10000 # 5000
+iObs = 2:4 # Indices for the 3 +ve observation types - eventually this will be replaced by an index for compartment I (and possibly C)
+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)
+radiusAS = c(1,2)
+(scenarios = expand.grid(flagFence=flagFence, flagHuntZone=flagHuntZone, radiusAS=radiusAS))
+nScenarios = nrow(scenarios)
+scenarioList = vector("list", nScenarios)
+### for (scenario in 1:nScenarios) { # scenario=1
+for (scenario in SCENARIO) { # scenario=1
+ nimPrint("################ Scenario=", scenario, " ################")
+ nimPrint(mcmcDir)
+ #####################
+ ## Simulate epidemics
+ #####################
+ mySims <- cSimFromMCMC$run(nSim, flagFence=scenarios[scenario,"flagFence"], flagHuntZone=scenarios[scenario,"flagHuntZone"], radiusAS=scenarios[scenario,"radiusAS"])
+ ##### Plot Maps #####
+ scenarioList[[scenario]] <- mySims
+ fileNames = c("Mean", "stDev", "MeanW", "stDevW", "Best")
+ for(day in seq(60, 140, by=20)) { # day=80
+ nimPrint(day)
+ for (iOut in 1:5) {
+ fileName = fileNames[iOut]
+ mapInfectiousDayX(mySims[iOut,,], Hex=hexAll,
+ Admin = admin %>%filter(is.element(rowAdmin, c(6,9,11,22))),
+ Day=day, outputDir=simDir,
+ fileName=paste0(fileName, "_sims",nSim, "_scenario", scenario,"_d",day),
+ StepsPerChunk=1, Fence=fence, HuntZone=huntZone,
+ pig_sites=pig_sites,
+ Outbreaks=outbreaks, caption=fileNames[iOut])
+ }
+ }
+ ##### Write output to files #####
+ (scenariosString = paste0(c("Fence", "Army","Radius"), scenarios[scenario,], collapse="_"))
+ for (iOut in 1:5) {
+ write.table(mySims[iOut,,],
+ col.names=FALSE,
+ row.names = paste0("D",(1:nrow(mySims[iOut,,]))-1),
+ file = paste0(simDir,"/",fileNames[iOut], "_", scenariosString, "_nSim", nSim, ".csv"))
+ }
+}
+
+
+
+##
+## (simFile = paste0(simDir, "/sims_nSim",nSim,"_scenario",SCENARIO,".Rdata"))
+## save(SCENARIO,
+## scenarios,
+## scenarioList,
+## admin,
+## bboxAll,
+## fence,
+## hexAll,
+## hexCentroidsAll,
+## outbreaks,
+## pig_sites,
+## stepsPerChunkAll,
+## wildBoarArrayAll,
+## wildBoarObsAll,
+## mapInfectiousDayX,
+## radiusActiveSearch,
+## file=simFile)
+
+
+## Write scenariosList as a series of files
+## load(simFile)
+
+
+## 02h58
+## Sys.time()
--
GitLab
From a7860446e10d4b90c4ba7238fb7efab2e047f9d6 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 11:44:08 +0100
Subject: [PATCH 3/9] Reverting model9 to it's form on the huntZone branch,
i.e. how we've used it.
---
nimble/model9/fitWildBoarModel9.R | 39 +++++++++++-----------------
nimble/model9/model_definition.R | 37 +++++++-------------------
nimble/model9/simWildBoarFromMCMC9.R | 39 ++++++++++------------------
src/plan.R | 7 +++--
4 files changed, 43 insertions(+), 79 deletions(-)
diff --git a/nimble/model9/fitWildBoarModel9.R b/nimble/model9/fitWildBoarModel9.R
index dcf581d..7526683 100644
--- a/nimble/model9/fitWildBoarModel9.R
+++ b/nimble/model9/fitWildBoarModel9.R
@@ -40,12 +40,8 @@ loadd(admin,
hexAll,
hexCentroidsSub,
fence,
- weightTauAHome1km,
- weightTauAAway1km,
- weightTauAHome2km,
- weightTauAAway2km,
- ## probNeighAS1km,
- ## probNeighAS2km,
+ probNeighAS1km,
+ probNeighAS2km,
outbreaks,
stepsPerChunkSub,
wildBoarArraySub, # observations aggregated over ten day intervals, hexagons restricted to 20km buffer
@@ -56,24 +52,19 @@ loadd(admin,
## Initial values and constants ####
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tStopSub <- 80
-initWildboarSub <- setWildBoarModelInitialValues(wbArray4sim = wildBoarArraySub,
- wildBoarObs = wildBoarObsSub,
- stepsPerChunk = stepsPerChunkSub,
- bboxAll = bboxAll) # For fitting MUST OMIT tStop ARGUMENT
-constWildboarSub <- setWildBoarModelConstants (wbArray4sim = wildBoarArraySub,
- Hex = hexSub,
- HexCentroids = hexCentroidsSub,
- Outbreaks = outbreaks,
- stepsPerChunk = stepsPerChunkSub,
- bboxAll = bboxAll,
- tStop = tStopSub,
- weightTauAHome1km = weightTauAHome1km,
- weightTauAAway1km = weightTauAAway1km,
- weightTauAHome2km = weightTauAHome2km,
- weightTauAAway2km = weightTauAAway2km)
- # probNeighAS1km = probNeighAS1km,
- # probNeighAS2km = probNeighAS2km
-
+initWildboarSub <- setWildBoarModelInitialValues(wbArray4sim = wildBoarArraySub,
+ wildBoarObs = wildBoarObsSub,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll) # For fitting MUST OMIT tStop ARGUMENT
+constWildboarSub <- setWildBoarModelConstants (wbArray4sim = wildBoarArraySub,
+ Hex = hexSub,
+ HexCentroids = hexCentroidsSub,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkSub,
+ bboxAll = bboxAll,
+ tStop = tStopSub,
+ probNeighAS1km = probNeighAS1km,
+ probNeighAS2km = probNeighAS2km)
##%%%%%%%%%%%%%%%%%%%
## Build R model ####
##%%%%%%%%%%%%%%%%%%%
diff --git a/nimble/model9/model_definition.R b/nimble/model9/model_definition.R
index 42d9a2e..31af076 100644
--- a/nimble/model9/model_definition.R
+++ b/nimble/model9/model_definition.R
@@ -58,9 +58,7 @@ bugsBoarCode = nimbleCode({
pRadius[1] <- 1.0
pRadius[2] <- 0.0
radiusActiveSearch ~ dcat(pRadius[1:2]) # Can set node value to either 1 or 2 in R code, to switch between two fixed parameters
- ## probNeighAS <- probNeighAS1km * (radiusActiveSearch==1) + probNeighAS2km * (radiusActiveSearch==2)
- weightTauAHome <- weightTauAHome1km * (radiusActiveSearch==1) + weightTauAHome2km * (radiusActiveSearch==2)
- weightTauAAway <- weightTauAAway1km * (radiusActiveSearch==1) + weightTauAAway2km * (radiusActiveSearch==2)
+ probNeighAS <- probNeighAS1km * (radiusActiveSearch==1) + probNeighAS2km * (radiusActiveSearch==2)
## Observed &/or Simulated Data
obsY[1:(nObsStates + returnI), 1:nTimeChunksPlus1, 1:nPops] <- simulateObservations(
tIntro = tIntro, xIntro = xIntro, yIntro = yIntro,
@@ -75,19 +73,14 @@ bugsBoarCode = nimbleCode({
neighStart = neighStart[1:nPops], neighEnd = neighEnd[1:nPops],
seperatedByFence = seperatedByFence[1:lNeighVec],
omegaFence = omegaFence, tFence = tFence, lNeighVec = lNeighVec,
- inHuntZone = inHuntZone[1:nPops],
- # probNeighAS = probNeighAS,
- weightTauAHome = weightTauAHome,
- weightTauAAway = weightTauAAway,
- radiusActiveSearch = radiusActiveSearch,
+ inHuntZone = inHuntZone[1:nPops], probNeighAS = probNeighAS,
tDelayActiveSearch = tDelayActiveSearch, decayTauA = decayTauA)
# For monitor2
loglik <- 0
})
-setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks, stepsPerChunk, bboxAll, tStop, returnI=FALSE,
- weightTauAHome1km, weightTauAAway1km, weightTauAHome2km, weightTauAAway2km) { ## probNeighAS1km, probNeighAS2km
+setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks, stepsPerChunk, bboxAll, tStop, returnI=FALSE, probNeighAS1km, probNeighAS2km) {
## wbArray4sim - an element of {wildBoarArrayAll, wildBoarArraySub}
xlims = bboxAll[c("xmin","xmax")]
ylims = bboxAll[c("ymin","ymax")]
@@ -124,12 +117,8 @@ setWildBoarModelConstants <- function (wbArray4sim, Hex, HexCentroids, Outbreaks
yMax = ylims[2],
seperatedByFence = attributes(Hex)$seperatedByFence,
inHuntZone = 1.0 * Hex$inHuntZone,
- weightTauAHome1km = weightTauAHome1km,
- weightTauAAway1km = weightTauAAway1km,
- weightTauAHome2km = weightTauAHome2km,
- weightTauAAway2km = weightTauAAway2km
- ## probNeighAS1km = probNeighAS1km,
- ## probNeighAS2km = probNeighAS2km
+ probNeighAS1km = probNeighAS1km,
+ probNeighAS2km = probNeighAS2km
## x0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"X"]), # X coordinate of first observed infected farm
## y0 = as.double((Outbreaks %>% filter(DATE.SUSP==0) %>% st_coordinates())[1,"Y"]) # Y coordinate of first observed infected farm
)
@@ -212,10 +201,7 @@ simulateObservations = nimbleFunction (
tFence = double(0, default=1E6),
lNeighVec = integer(0),
inHuntZone = double(1),
- ## probNeighAS = double(0),
- weightTauAHome = double(0),
- weightTauAAway = double(0),
- radiusActiveSearch = integer(0),
+ probNeighAS = double(0),
tDelayActiveSearch = integer(0),
decayTauA = double(0) ) {
# Indices for states
@@ -310,12 +296,7 @@ simulateObservations = nimbleFunction (
pHome_x_pAway_x_2 = 2*pHome*pAway
pAway_x_pAway_x_2 = 2*pAway*pAway
nNeigh_x_nNeigh = nimNumeric(length=nPops, value=nNeigh[1:nPops]*nNeigh[1:nPops])
- ## Adjust tauA for home/away pixels and 1km/2km radius
- ## i.e. the weighting affects the daily probability, not the rate.
- tauAHome = -log(1 - weightTauAHome * (1 - exp(-tauA)))
- tauAAway = -log(1 - weightTauAAway * (1 - exp(-tauA)))
- ## probNeighAS_x_tauA = probNeighAS * tauA
- ##
+ probNeighAS_x_tauA = probNeighAS * tauA
## Dynamics
for (t in floor(tIntro):tStop) {
# Indicator to turn off/on observation processes
@@ -420,10 +401,10 @@ simulateObservations = nimbleFunction (
if (indicatorUpdateAS & (dStates[Hpos] + dStates[Sa] + dStates[Sp] > 0)) {
# Then an active search will start in 3 days (delay for test) + 7 days (fixed delay for organising search in simulated data)
tauA_matrix[pop, tPlus1] = decayTauA * tauA_matrix[pop, t]
- tauA_matrix[pop, tPlusDelay] = tauAHome
+ tauA_matrix[pop, tPlusDelay] = tauA
for (ii in neighStart[pop]:neighEnd[pop]) {
iNeigh = neighVec[ii]
- tauA_matrix[iNeigh, tPlusDelay] = min(tauAHome, tauAAway + tauA_matrix[iNeigh, tPlusDelay])
+ tauA_matrix[iNeigh, tPlusDelay] = min(tauA, probNeighAS_x_tauA + tauA_matrix[iNeigh, tPlusDelay])
}
}
} # End loop on pop
diff --git a/nimble/model9/simWildBoarFromMCMC9.R b/nimble/model9/simWildBoarFromMCMC9.R
index 1c9568e..087628a 100644
--- a/nimble/model9/simWildBoarFromMCMC9.R
+++ b/nimble/model9/simWildBoarFromMCMC9.R
@@ -1,7 +1,7 @@
# source(here::here("nimble/model9/simWildBoarFromMCMC9.R"))
# rm(list=ls())
-SCENARIO <- 1
+SCENARIO <- 8
# This script jumps through MCMC output and generates weighted means & stdevs of the simulations
@@ -43,12 +43,8 @@ if (reBuildModels) {
huntZone,
hexAll,
hexCentroidsAll,
- weightTauAHome1km,
- weightTauAAway1km,
- weightTauAHome2km,
- weightTauAAway2km,
- #probNeighAS1km,
- #probNeighAS2km,
+ probNeighAS1km,
+ probNeighAS2km,
outbreaks,
pig_sites,
stepsPerChunkAll,
@@ -64,20 +60,16 @@ if (reBuildModels) {
bboxAll = bboxAll,
tStop = tStopAll,
returnI = TRUE)
- constWildboarAll = setWildBoarModelConstants(wbArray4sim = wildBoarArrayAll,
- Hex = hexAll,
- HexCentroids = hexCentroidsAll,
- Outbreaks = outbreaks,
- stepsPerChunk = stepsPerChunkAll,
- bboxAll = bboxAll,
- tStop = tStopAll,
- weightTauAHome1km = weightTauAHome1km,
- weightTauAAway1km = weightTauAAway1km,
- weightTauAHome2km = weightTauAHome2km,
- weightTauAAway2km = weightTauAAway2km,
- ## probNeighAS1km = probNeighAS1km,
- ## probNeighAS2km = probNeighAS2km,
- returnI = TRUE)
+ constWildboarAll = setWildBoarModelConstants(wbArray4sim = wildBoarArrayAll,
+ Hex = hexAll,
+ HexCentroids = hexCentroidsAll,
+ Outbreaks = outbreaks,
+ stepsPerChunk = stepsPerChunkAll,
+ bboxAll = bboxAll,
+ tStop = tStopAll,
+ probNeighAS1km = probNeighAS1km,
+ probNeighAS2km = probNeighAS2km,
+ returnI = TRUE)
##%%%%%%%%%%%%%%%%%%%
## Build R model ####
##%%%%%%%%%%%%%%%%%%%
@@ -224,7 +216,4 @@ for (scenario in SCENARIO) { # scenario=1
## Write scenariosList as a series of files
## load(simFile)
-
-
-## 02h58
-## Sys.time()
+## probNeighAS = probNeighAS1km
diff --git a/src/plan.R b/src/plan.R
index 54b90d4..04e9415 100644
--- a/src/plan.R
+++ b/src/plan.R
@@ -598,12 +598,15 @@ plan <- drake_plan(
innerOuterAreas2km = sampleProbActive(hexTiny, cellSize, radiusInM=2000, resolution=100)
,
+ ## New weighting for model10
weightTauAHome1km = 1, ## innerOuterAreas1km[1]/innerOuterAreas1km[1]
weightTauAAway1km = innerOuterAreas1km[2] / innerOuterAreas1km[1],
weightTauAHome2km = innerOuterAreas2km[1] / innerOuterAreas1km[1],
weightTauAAway2km = innerOuterAreas2km[2] / innerOuterAreas1km[1],
- ## probNeighAS2km =
- ## probNeighAS1km =
+
+ ## Old weighting for model9 (and previous).
+ probNeighAS1km = weightTauAAway1km,
+ probNeighAS2km = weightTauAAway2km,
## hexWhat = st_make_grid(hexTiny, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
## hexWhat = st_make_grid(HexCenter, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
--
GitLab
From 9ab5319380c25d5d4937678af20812cc6e0c2456 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 11:55:44 +0100
Subject: [PATCH 4/9] Importing some changes on branch improved_active_search.
---
.gitignore | 2 --
src/functions.R | 12 ++++++++----
src/plan.R | 14 ++++++++++++--
3 files changed, 20 insertions(+), 8 deletions(-)
diff --git a/.gitignore b/.gitignore
index cfe85aa..49c5781 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,8 +20,6 @@ List of transmission pathways and risks.docx
gadm36_FRA_2_sp.rds
biblio
mostRecentMCMCSamples_BU*.csv
-MCMC
-MCMC1
wildBoarFitEpidemic_APT.R
asf-migale.tar
*Rdata
diff --git a/src/functions.R b/src/functions.R
index 208e9c1..80ede0b 100644
--- a/src/functions.R
+++ b/src/functions.R
@@ -1404,7 +1404,7 @@ makeHexTiny <- function(HexAll, CellSize) { # HexAll=hexAll; CellSize=cellSize
return(HexTiny)
}
-sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) { # HexTiny=hexTiny; CellSize=cellSize
+sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=100) { # HexTiny=hexTiny; CellSize=cellSize
radiusM <- set_units(radiusInM, m)
centroids <- st_centroid(HexTiny) %>% st_geometry()
distMat <- st_distance(centroids, centroids)
@@ -1412,9 +1412,10 @@ sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) {
HexCenter <- HexTiny[iCenter,] %>% st_geometry()
HexOuter <- st_difference(HexTiny, HexCenter) %>% select(geometry)
##
- delta <- CellSize * 0.6
+ delta <- CellSize * 0.6
##
grid = st_make_grid(HexCenter, cellsize = resolution, square = FALSE, offset = c(0, 0))
+ ## plot(HexTiny[,1])
## plot(grid, add=TRUE, col="blue")
dGrid <- data.frame(id=1:length(grid)) # data.frame(id=1:length(grid))
dGrid$geometry = grid # grid
@@ -1434,8 +1435,11 @@ sampleProbActive <- function(HexTiny, CellSize, radiusInM=1000, resolution=10) {
##
sumOuter <- overlapsOuter %>% st_area() %>% sum()
sumInner <- overlapsInner %>% st_area() %>% sum()
- (proportionActiveSearchInOneNeighbour <- sumOuter / sumInner)
- return(proportionActiveSearchInOneNeighbour)
+ ## Original output
+ ## (proportionActiveSearchInOneNeighbour <- sumOuter / sumInner)
+ ## return(proportionActiveSearchInOneNeighbour)
+ ## New output
+ return(c(sumInner, sumOuter))
}
diff --git a/src/plan.R b/src/plan.R
index 1765609..04e9415 100644
--- a/src/plan.R
+++ b/src/plan.R
@@ -593,11 +593,21 @@ plan <- drake_plan(
## They integrate over all possible locations (within the central hexagon) for an Active Search centroid AND over all possible locations visited by the active search
## This provides a weight for the strength of an active search in a pixel that does not containing the active search centroid (i.e. a neighbour)
## These weights will be added as constants in model8 to reduce bias from the current (false) assumption that the active search is globally homogeneous
- probNeighAS1km = sampleProbActive(hexTiny, cellSize, radiusInM=1000, resolution=100)
+ innerOuterAreas1km = sampleProbActive(hexTiny, cellSize, radiusInM=1000, resolution=100)
,
- probNeighAS2km = sampleProbActive(hexTiny, cellSize, radiusInM=2000, resolution=100)
+ innerOuterAreas2km = sampleProbActive(hexTiny, cellSize, radiusInM=2000, resolution=100)
,
+ ## New weighting for model10
+ weightTauAHome1km = 1, ## innerOuterAreas1km[1]/innerOuterAreas1km[1]
+ weightTauAAway1km = innerOuterAreas1km[2] / innerOuterAreas1km[1],
+ weightTauAHome2km = innerOuterAreas2km[1] / innerOuterAreas1km[1],
+ weightTauAAway2km = innerOuterAreas2km[2] / innerOuterAreas1km[1],
+
+ ## Old weighting for model9 (and previous).
+ probNeighAS1km = weightTauAAway1km,
+ probNeighAS2km = weightTauAAway2km,
+
## hexWhat = st_make_grid(hexTiny, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
## hexWhat = st_make_grid(HexCenter, n=1, cellsize = cellSize + 1000, square = FALSE, offset = c(0, 0))
## hexWhat <- hexWhat[1,]
--
GitLab
From f0c9c7a35d14f6da7656f3e0aceb0605f2ae742f Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 12:04:12 +0100
Subject: [PATCH 5/9] Updating file path names in model10 (i.e. changing 9s to
10s).
---
...itWildBoarModel9.R => fitWildBoarModel10.R} | 18 +++++++++---------
nimble/model10/model_definition.R | 5 -----
...BoarFromMCMC9.R => simWildBoarFromMCMC10.R} | 8 ++++----
3 files changed, 13 insertions(+), 18 deletions(-)
rename nimble/model10/{fitWildBoarModel9.R => fitWildBoarModel10.R} (97%)
rename nimble/model10/{simWildBoarFromMCMC9.R => simWildBoarFromMCMC10.R} (95%)
diff --git a/nimble/model10/fitWildBoarModel9.R b/nimble/model10/fitWildBoarModel10.R
similarity index 97%
rename from nimble/model10/fitWildBoarModel9.R
rename to nimble/model10/fitWildBoarModel10.R
index dcf581d..7842bd4 100644
--- a/nimble/model10/fitWildBoarModel9.R
+++ b/nimble/model10/fitWildBoarModel10.R
@@ -1,11 +1,11 @@
-# source(here::here("nimble/model9/fitWildBoarModel9.R"))
+# source(here::here("nimble/model10/fitWildBoarModel10.R"))
# rm(list=ls())
-iNSim = 4
+iNSim = 1
source(here::here("src/packages.R"))
source(here::here("src/functions.R"))
source(here::here("nimble/nimbleFunctions.R"))
-source(here::here("nimble/model9/model_definition.R"))
+source(here::here("nimble/model10/model_definition.R"))
options(width=1000)
# drake::r_make() # TO UPDATE INIT AND CONST
@@ -16,12 +16,12 @@ options(width=1000)
## nimPrint("iAttractI=",iAttractI, " logit(attractI)=",logit_attractI)
-(baseDir <- here::here("nimble/model9"))
+(baseDir <- here::here("nimble/model10"))
setwd(baseDir)
## (mcmcDir <- paste0(baseDir, "/MCMC_lpAI",logit_attractI))
-nSimVec = c(10, 50, 100, 500, 1000, 2000)
+nSimVec = c(500)
nSim = nSimVec[iNSim]
(mcmcDir = paste0(baseDir, paste0("/MCMC",1:length(nSimVec))[iNSim])) ## "/MCMC2"))
@@ -199,12 +199,12 @@ if (TRUE) { # FALSE
nimPrint(tenDayChunk, " ", iCol)
mapInfectiousDayX(t(rWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
- Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model10/checks",
fileName=paste0("check_chunk", tenDayChunk, "_", fNames[iCol], "_data"),
StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
mapInfectiousDayX(t(cWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
- Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model10/checks",
fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_sim"),
StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
}
@@ -339,12 +339,12 @@ for (iRun in 1:100) { # iRun=1
nimPrint(tenDayChunk, " ", iCol)
mapInfectiousDayX(t(wildBoarArraySub[iCol,,]), Hex=hexSub,
Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
- Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model10/checks",
fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_data"),
StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
mapInfectiousDayX(t(cWildBoarModelSub$obsY[iCol,,]), Hex=hexSub,
Admin=admin%>%filter(is.element(rowAdmin, c(6,9,11,22))),
- Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model9/checks",
+ Day=tenDayChunk, outputDir="/home/david/asf-challenge/nimble/model10/checks",
fileName=paste0("chunk", tenDayChunk, "_", fNames[iCol], "_sim"),
StepsPerChunk=10, Fence=fence, pig_sites=pig_sites, Outbreaks=outbreaks)
}
diff --git a/nimble/model10/model_definition.R b/nimble/model10/model_definition.R
index 42d9a2e..c0cd628 100644
--- a/nimble/model10/model_definition.R
+++ b/nimble/model10/model_definition.R
@@ -1,8 +1,3 @@
-## TO-DO
-## - visualise relation between passive and active search
-## - adapt to hexagonal grid
-
-
bugsCode4MV = nimbleCode({
logit(scTIntro) ~ dLogitBeta(7, 7) # Introduction date (on 0:1)
logit(pIntroX) ~ dLogitBeta(1, 1) # Introduction location relative to bounding box of whole island
diff --git a/nimble/model10/simWildBoarFromMCMC9.R b/nimble/model10/simWildBoarFromMCMC10.R
similarity index 95%
rename from nimble/model10/simWildBoarFromMCMC9.R
rename to nimble/model10/simWildBoarFromMCMC10.R
index 1c9568e..0fe181f 100644
--- a/nimble/model10/simWildBoarFromMCMC9.R
+++ b/nimble/model10/simWildBoarFromMCMC10.R
@@ -1,4 +1,4 @@
-# source(here::here("nimble/model9/simWildBoarFromMCMC9.R"))
+# source(here::here("nimble/model10/simWildBoarFromMCMC10.R"))
# rm(list=ls())
SCENARIO <- 1
@@ -10,9 +10,9 @@ SCENARIO <- 1
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
library(here)
## help.start()
-(modelDir = here::here("nimble/model9")) # For selecting model version
-(mcmcDir = paste0(modelDir,"/MCMC4")) # 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,"/sims")) # For selecting model version
+(modelDir = here::here("nimble/model10")) # For selecting model version
+(mcmcDir = paste0(modelDir,"/MCMC1")) # 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,"/sims")) # For selecting model version
setwd(modelDir)
# setwd(mcmcDir)
mcmcFile = paste0(mcmcDir,"/samples.csv") ## Assumed to be parameters
--
GitLab
From e6faeb63596bddf14a42f295fda79e54f53d33b8 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 12:13:32 +0100
Subject: [PATCH 6/9] .gitignore update
---
.gitignore | 2 --
1 file changed, 2 deletions(-)
diff --git a/.gitignore b/.gitignore
index cfe85aa..49c5781 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,8 +20,6 @@ List of transmission pathways and risks.docx
gadm36_FRA_2_sp.rds
biblio
mostRecentMCMCSamples_BU*.csv
-MCMC
-MCMC1
wildBoarFitEpidemic_APT.R
asf-migale.tar
*Rdata
--
GitLab
From 749743a9576eecd57b435f150f3daf974d9b39cb Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 12:23:20 +0100
Subject: [PATCH 7/9] commit-size pearl script - shows the disk space usage of
each commit.
---
.gitignore | 2 ++
src/commit-size | 16 ++++++++++++++++
2 files changed, 18 insertions(+)
create mode 100755 src/commit-size
diff --git a/.gitignore b/.gitignore
index 49c5781..24abd48 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,3 +31,5 @@ sims
tar_*
*.tar
simple-feature-example.R
+model7
+model8
diff --git a/src/commit-size b/src/commit-size
new file mode 100755
index 0000000..91323ec
--- /dev/null
+++ b/src/commit-size
@@ -0,0 +1,16 @@
+#!/usr/bin/perl
+foreach my $rev (`git rev-list --all --pretty=oneline`) {
+ my $tot = 0;
+ ($sha = $rev) =~ s/\s.*$//;
+ foreach my $blob (`git diff-tree -r -c -M -C --no-commit-id $sha`) {
+ $blob = (split /\s/, $blob)[3];
+ next if $blob == "0000000000000000000000000000000000000000"; # Deleted
+ my $size = `echo $blob | git cat-file --batch-check`;
+ $size = (split /\s/, $size)[2];
+ $tot += int($size);
+ }
+ my $revn = substr($rev, 0, 40);
+# if ($tot > 1000000) {
+ print "$tot $revn " . `git show --pretty="format:" --name-only $revn | wc -l` ;
+# }
+}
--
GitLab
From 9ce36a3a82070546f0e01efbf71b91f493aa75a2 Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 12:50:49 +0100
Subject: [PATCH 8/9] Edits in fitWildBoarModel10.R
---
.gitignore | 8 ++++++++
nimble/model10/fitWildBoarModel10.R | 11 ++++++-----
2 files changed, 14 insertions(+), 5 deletions(-)
diff --git a/.gitignore b/.gitignore
index 24abd48..03e9660 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,11 @@ tar_*
simple-feature-example.R
model7
model8
+BU
+BU-sunday
+BU-sunday-night
+model9/MCMC1/*csv
+model9/MCMC2/*csv
+model9/MCMC5/*csv
+model9/MCMC6/*csv
+model9/MCMC7/*csv
diff --git a/nimble/model10/fitWildBoarModel10.R b/nimble/model10/fitWildBoarModel10.R
index 7842bd4..8928b7d 100644
--- a/nimble/model10/fitWildBoarModel10.R
+++ b/nimble/model10/fitWildBoarModel10.R
@@ -1,6 +1,6 @@
# source(here::here("nimble/model10/fitWildBoarModel10.R"))
-# rm(list=ls())
-iNSim = 1
+rm(list=ls()) # If sourcing this will ensure the cModel is rebuilt (needed because line 172 could be dangerous)
+iNSim = 5
source(here::here("src/packages.R"))
source(here::here("src/functions.R"))
@@ -21,7 +21,7 @@ setwd(baseDir)
## (mcmcDir <- paste0(baseDir, "/MCMC_lpAI",logit_attractI))
-nSimVec = c(500)
+nSimVec = c(500, 500, 500, 500, 500)
nSim = nSimVec[iNSim]
(mcmcDir = paste0(baseDir, paste0("/MCMC",1:length(nSimVec))[iNSim])) ## "/MCMC2"))
@@ -98,7 +98,7 @@ stochNodes = rWildBoarModelSub$getNodeNames(stochOnly=TRUE)
#################################
Fixed <- c("logit_pR", "radiusActiveSearch")
Target <- sort(c(
- "logit_pHuntY",
+ "logit_pHuntY", ## Not good idea to fix this - it is too low.
"logit_scTIntro",
"logit_pIntroX",
"logit_pIntroY",
@@ -146,7 +146,6 @@ if (TRUE) { # FALSE
(rWildBoarModelSub$log_beta <- previousSamples$log_beta)
(iniPropCov <- diag(1E-3, length(Target)))
}
- (rWildBoarModelSub$logit_pHuntY <- 0) # previousSamples$logit_pHuntY)
(rWildBoarModelSub$logit_pR <- logit(0.05))
}
@@ -170,8 +169,10 @@ for (tt in Target) {
## Compile models ####
##%%%%%%%%%%%%%%%%%%%%
if (is.element("cWildBoarModelSub", ls())) {
+ nimPrint("WARNING: NOT RECOMPILING cWildBoarModelSub. USING PARAMETERS FROM rWildBoarModelSub. VERIFY THAT THIS IS DESIRED BEHAVIOUR!!!")
nimCopy(from=rWildBoarModelSub, to=cWildBoarModelSub, nodes = stochNodes)
} else {
+ nimPrint("RECOMPILING cWildBoarModelSub")
cWildBoarModelSub = compileNimble(rWildBoarModelSub, showCompilerOutput = TRUE) # FALSE
}
--
GitLab
From 0a4fbc75c571b67e59241e6f76bfaa74620b104c Mon Sep 17 00:00:00 2001
From: david <david.pleydell@inrae.fr>
Date: Mon, 23 Nov 2020 13:15:26 +0100
Subject: [PATCH 9/9] Updated paths for setting initial values in model10
MCMCs.
---
nimble/model10/fitWildBoarModel10.R | 14 +++++++-------
1 file changed, 7 insertions(+), 7 deletions(-)
diff --git a/nimble/model10/fitWildBoarModel10.R b/nimble/model10/fitWildBoarModel10.R
index 8928b7d..91c187d 100644
--- a/nimble/model10/fitWildBoarModel10.R
+++ b/nimble/model10/fitWildBoarModel10.R
@@ -129,20 +129,20 @@ if (TRUE) { # FALSE
nimPrint(Target[iTarget], " = ", rWildBoarModelSub[[Target[iTarget]]])
}
} else {
- model1SamplesFile <- "/home/david/asf-challenge/nimble/model1/MCMC/samples.csv"
+ model1SamplesFile <- "/home/david/asf-challenge/nimble/model9/MCMC4/samples.csv"
nimPrint("Using ", model1SamplesFile)
previousSamples <- read.table(file=model1SamplesFile, header = TRUE)
previousSamples <- tail(previousSamples,1)
- (rWildBoarModelSub$logit_scTIntro <- logit(1-(14)/100))
+ (rWildBoarModelSub$logit_scTIntro <- previousSamples$logit_scTIntro)
(rWildBoarModelSub$logit_pIntroX <- previousSamples$logit_pIntroX)
(rWildBoarModelSub$logit_pIntroY <- previousSamples$logit_pIntroY)
(rWildBoarModelSub$logit_pHome <- previousSamples$logit_pHome)
(rWildBoarModelSub$logit_attractI <- previousSamples$logit_attractI)
- (rWildBoarModelSub$logit_omegaFence <- logit(0.999))
- (rWildBoarModelSub$log_tauP <- -8) # 1 + log(exp(previousSamples$log_tauCarDet) * (1-ilogit(previousSamples$logit_pActive))))
- (rWildBoarModelSub$log_tauPArmy <- -6.6)
- (rWildBoarModelSub$log_tauA <- 2 + log(exp(previousSamples$log_tauCarDet) * ilogit(previousSamples$logit_pActive)))
- (rWildBoarModelSub$log_tauHArmy <- -3.47)
+ (rWildBoarModelSub$logit_omegaFence <- previousSamples$logit_omegaFence)
+ (rWildBoarModelSub$log_tauP <- previousSamples$log_tauP)
+ (rWildBoarModelSub$log_tauPArmy <- previousSamples$log_tauPArmy)
+ (rWildBoarModelSub$log_tauA <- previousSamples$log_tauA)
+ (rWildBoarModelSub$log_tauHArmy <- previousSamples$log_tauHArmy)
(rWildBoarModelSub$log_beta <- previousSamples$log_beta)
(iniPropCov <- diag(1E-3, length(Target)))
}
--
GitLab