########################################
#
#   MoNAn workshop Oct 2023, Zurich
#   Script 2
#
########################################

### this script goes through some potential reasons
# why a model might not achieve convergence
# (see the manual for further information)

library(MoNAn)

##### data preparation based on script 1 #####

N_ind <- nrow(mobilityEdgelist)
N_org <- length(unique(as.numeric(mobilityEdgelist)))
transfers <- createEdgelist(mobilityEdgelist,
                            nodeSet = c("organisations", "organisations", "people")
)
people <- createNodeSet(N_ind)
organisations <- createNodeSet(N_org)
region <- createNodeVariable(orgRegion, nodeSet = "organisations", addSame = T)
sex <- createNodeVariable(indSex, nodeSet = "people")
# combine created objects to the process state
myState <- createProcessState(
  list(
    transfers = transfers,
    people = people,
    organisations = organisations,
    region = region,
    sex = sex), 
  dependentVariable = "transfers")
myCache <- createWeightedCache(myState, resourceCovariates = c("sex"))



##### Degenerate model #####

# effects object
myEffects_degen <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_basic"),
    list("same_covariate", attribute.index = "region"),
    list("loops_resource_covar", resource.attribute.index = "sex")
  ), 
  checkProcessState = myState  
)

myAlg_degen <- createAlgorithm(myState, myEffects_degen, multinomialProposal = FALSE,
                               initGain = 0.8,
                               initialIterationsN2 = 70)

myResDN_degen <- estimateMobilityNetwork(
  myState, myCache, myEffects_degen, myAlg_degen,
  parallel = TRUE, cpus = 4,
  verbose = T,
  returnDeps = T
)

myResDN_degen

myAlg_degen2 <- createAlgorithm(myState, myEffects_degen, multinomialProposal = FALSE,
                                initGain = 0,
                                nsubN2 = 1,
                                initialIterationsN2 = 1,
                                iterationsN3 = 1000)

myResDN_degen2 <- estimateMobilityNetwork(
  myState, myCache, myEffects_degen, myAlg_degen2,
  initialParameters = myResDN_degen$estimates,
  parallel = TRUE, cpus = 2,
  verbose = T,
  returnDeps = T
)

traces_degen <- extractTraces(myResDN_degen2, myEffects_degen)

plot(traces_degen)


##### for non-degenerate models #####

##### the baseline (converged) specification #####

# effects object
myEffects <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_min"),
    list("same_covariate", attribute.index = "region"),
    list("loops_resource_covar", resource.attribute.index = "sex")
  ), 
  checkProcessState = myState  
)

myAlg <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                         initGain = 0.8)

myResDN <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg,
  parallel = TRUE, cpus = 4,
  verbose = TRUE
)

myResDN

##### Too small thinning in Phase 2 #####

myAlg_thin_N2 <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                           thinningN2 = 20)

myResDN_thin_N2 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_thin_N2,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_thin_N2

# solution: check for thinning in phase 3 and adjust to the checked value

##### Too small thinning in Phase 3 #####

myAlg_thin_N3 <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                           nsubN2 = 1,
                           initialIterationsN2 = 1,
                           initGain = 0,
                           thinningN3 = 100)

myResDN_thin_N3 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_thin_N3,
  initialParameters = myResDN$estimates,
  parallel = TRUE, cpus = 4,
  verbose = T,
  returnDeps = T
)

myResDN_thin_N3

traces_thin_N3 <- extractTraces(myResDN_thin_N3, myEffects)

plot(traces_thin_N3)

# the plot shows auto-correlation, which indicates the thinning is too low


##### Too few steps in Phase 2 #####

myAlg_few_steps_N2 <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                           initialIterationsN2 = 10)

myResDN_few_steps_N2 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_few_steps_N2,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_few_steps_N2

myResDN_few_steps_N2_2 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_few_steps_N2,
  prevAns = myResDN_few_steps_N2,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_few_steps_N2
myResDN_few_steps_N2_2

# potentially the case if parameters change between subsequent runs


##### Too low gain in Phase 2 #####

myAlg_d <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                           initGain = 0.02)

myResDN_d_1 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_d,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_d_2 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_d,
  prevAns = myResDN_d_1,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_d_1
myResDN_d_2

# potentially the case if parameters change between subsequent runs


##### Too few steps in Phase 3 #####

myAlg_few_steps_N3 <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                                      nsubN2 = 1,
                                      initialIterationsN2 = 1,
                                      initGain = 0,
                                      iterationsN3 = 100)

myResDN_few_steps_N3_1 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_few_steps_N3,
  initialParameters = myResDN$estimates,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_few_steps_N3_2 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_few_steps_N3,
  initialParameters = myResDN$estimates,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_few_steps_N3_3 <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_few_steps_N3,
  initialParameters = myResDN$estimates,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_few_steps_N3_1
myResDN_few_steps_N3_2
myResDN_few_steps_N3_3

# parameters remain but convergence stats change


##### BONUS #####

##### Too few steps in Phase 1 #####

myAlg_a <- createAlgorithm(myState, myEffects, multinomialProposal = FALSE,
                         burnInN1 = 100,
                         thinningN1 = 50,
                         iterationsN1 = 3)

myResDN_a <- estimateMobilityNetwork(
  myState, myCache, myEffects, myAlg_a,
  parallel = TRUE, cpus = 4,
  verbose = T
)

myResDN_a

# the simplest solution here is not to deviate from the default values for
# Phase 1





