## ----setup, message = FALSE, warning = FALSE----------------------------------
library(sbm)
library(ggplot2)
library(patchwork) # to put figures side by side
theme_set(theme_bw())


## ----show function------------------------------------------------------------
showSBM <- function(mySBM, title) {
  (plot(mySBM, type = "data") + ggtitle(title)) |
    plot(mySBM, type = "graphon") |
    plot(mySBM, type = "meso", plotOptions = list(seed = 1))
}


## ----communities--------------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
set.seed(1)
alphaCommunities <- matrix(0.05, 3, 3) + diag(0.5, 3)
alphaCommunities
communities <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.4, 0.35, 0.25),
                               connectParam = list(mean = alphaCommunities))
showSBM(communities, "Communities")


## ----core-periphery-----------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaCore <- matrix(c(0.9, 0.6, 0.3,
                      0.6, 0.2, 0.05,
                      0.3, 0.05, 0.01), 3, 3)
corePeriphery <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.15, 0.35, 0.5),
                                 connectParam = list(mean = alphaCore))
showSBM(corePeriphery, "Core and periphery")


## ----disassortative-----------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaDisassortative <- matrix(c(0.05, 0.6,
                                0.6, 0.05), 2, 2)
disassortative <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.5, 0.5),
                                  connectParam = list(mean = alphaDisassortative))
showSBM(disassortative, "Disassortative")


## ----hierarchy----------------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaHierarchy <- matrix(c(0.2, 0.6, 0.3,
                           0.02, 0.2, 0.7,
                           0.01, 0.02, 0.2), 3, 3, byrow = TRUE)
hierarchy <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.3, 0.3, 0.4),
                             connectParam = list(mean = alphaHierarchy),
                             directed = TRUE)
showSBM(hierarchy, "Directed hierarchy")


## ----poisson------------------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaPoisson <- matrix(c(8, 2, 0.5,
                         2, 4, 1,
                         0.5, 1, 0.2), 3, 3)
countNetwork <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.3, 0.3, 0.4),
                                connectParam = list(mean = alphaPoisson),
                                model = "poisson")
showSBM(countNetwork, "Poisson SBM")


## ----gaussian-----------------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaGaussian <- matrix(c(4, 1, 0,
                          1, 3, 2,
                          0, 2, 1), 3, 3)
weightedNetwork <- sampleSimpleSBM(nbNodes = 90, blockProp = c(0.3, 0.3, 0.4),
                                   connectParam = list(mean = alphaGaussian, var = 1),
                                   model = "gaussian")
showSBM(weightedNetwork, "Gaussian SBM")


## ----bipartite modular--------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaModular <- matrix(c(0.7, 0.05, 0.05,
                         0.05, 0.6, 0.05,
                         0.05, 0.05, 0.8), 3, 3)
modular <- sampleBipartiteSBM(nbNodes = c(60, 90),
                              blockProp = list(row = c(0.3, 0.4, 0.3), col = c(0.3, 0.3, 0.4)),
                              connectParam = list(mean = alphaModular),
                              dimLabels = c(row = "plant", col = "pollinator"))
showSBM(modular, "Modular bipartite network")


## ----bipartite nested---------------------------------------------------------
#| fig-width: 12
#| fig-height: 4
alphaNested <- matrix(c(0.95, 0.8, 0.5,
                        0.8, 0.4, 0.1,
                        0.5, 0.1, 0.02), 3, 3)
nested <- sampleBipartiteSBM(nbNodes = c(60, 90),
                             blockProp = list(row = c(0.2, 0.3, 0.5), col = c(0.2, 0.3, 0.5)),
                             connectParam = list(mean = alphaNested),
                             dimLabels = c(row = "plant", col = "pollinator"))
showSBM(nested, "Nested bipartite network")


## ----fit----------------------------------------------------------------------
fitCore <- estimateSimpleSBM(corePeriphery$networkData, model = "bernoulli",
                             estimOptions = list(verbosity = 0, plot = FALSE))


## ----ICL----------------------------------------------------------------------
fitCore$storedModels |>
  ggplot() + aes(x = nbBlocks, y = ICL) + geom_line() + geom_point() +
  geom_vline(xintercept = fitCore$nbBlocks, color = "red", linetype = "dashed")


## ----compare blocks-----------------------------------------------------------
table(estimated = fitCore$memberships, true = corePeriphery$memberships)
aricode::ARI(fitCore$memberships, corePeriphery$memberships)


## ----estimated parameters-----------------------------------------------------
round(fitCore$connectParam$mean, 2)
round(fitCore$blockProp, 2)


## ----compare meso-------------------------------------------------------------
#| fig-width: 10
triangle <- cbind(x = c(0, -1, 1), y = c(1, -0.7, -0.7))
plot(corePeriphery, type = "meso", plotOptions = list(layout = triangle, title = "True")) |
  plot(fitCore, type = "meso", plotOptions = list(layout = triangle, title = "Estimated"))


## ----data expected------------------------------------------------------------
#| fig-width: 10
plot(fitCore, type = "data") | plot(fitCore, type = "expected")


## ----meso options-------------------------------------------------------------
#| fig-width: 10
plot(fitCore, type = "meso", plotOptions = list(seed = 1, withinBlock = "fill")) |
  plot(fitCore, type = "meso", plotOptions = list(seed = 1, edge.threshold = 0.1,
                                                  vertex.color = "steelblue"))