##
## Module 11 helper: the cold-start bootstrap for the dependent multi-layer model.
##
## You source this file once in Part 2 of the tutorial and call init_cross_degree().
## You do not need to edit it; read it to see what the bootstrap does, then get back
## to the epidemiology. The per-step callback that keeps the cross-layer degree
## current during the simulation is short enough that the tutorial shows it inline.
##
## The cross-layer bootstrap below originates with Chad Klumb (EpiModel Gallery
## "multinets" example, co-authored with Samuel Jenness).
##

## init_cross_degree(): solve the cold-start (chicken-and-egg) problem.
##
## Each layer's formation model wants a "does this node have a partner in the OTHER
## layer" attribute (deg1+.net1 / deg1+.net2). But before we fit anything, neither
## layer exists yet, so those attributes have no values. We break the deadlock with
## san() (simulated annealing from ergm): draw a plausible edge set for each layer
## that hits its target statistics, read the degrees off it, and store them as the
## starting cross-layer attributes. These are only starting values; netest refines
## the models from here.
##
## Targets below match Part 2 of the tutorial (layer 1: 90 edges, 60 same-race;
## layer 2: 75 edges, 120 degree-1 nodes, 10 cross-layer). Returns the network with
## deg1+.net1 and deg1+.net2 set.
init_cross_degree <- function(nw) {
  # Draw layer 1 (main) and record who has a main partner.
  nw_san <- san(nw ~ edges + nodematch("race"), target.stats = c(90, 60))
  nw_san <- set_vertex_attribute(nw_san, "deg1+.net1",
                                 pmin(get_degree(nw_san), 1))

  # Draw layer 2 (casual) conditional on layer 1: the nodefactor term makes nodes
  # already active in layer 1 less likely to pick up a casual partner.
  nw_san_2 <- san(nw_san ~ edges + degree(1) + nodefactor("deg1+.net1"),
                  target.stats = c(75, 120, 10))

  # Store both starting cross-layer attributes on the network netest will use.
  nw <- set_vertex_attribute(nw, "deg1+.net1", pmin(get_degree(nw_san),   1))
  nw <- set_vertex_attribute(nw, "deg1+.net2", pmin(get_degree(nw_san_2), 1))
  nw
}