
# ABOUT -------------------------------------------------------------------

# This is a replication script of the paper blockmodeling: an R package for
# Generalized Blockmodeling, published in the Journal of Statistical 
# Software.
# 
# Authors: Miha Matjašič (miha.matjasic@fdv.uni-lj.si), Marjan Cugmas, Aleš Žiberna.


# SECTION 4.: Demonstration of the package use ----------------------------

# Package
library(blockmodeling)

# Data
data("baker")

# reading from  text file
baker<-read.delim("baker original selfties.txt",row.names = 1)
#reading from Pajek .net file
baker<-loadnetwork("baker_selfties.net")

#remove diagonal
diag(baker) <- 0

# Visualization
plotMat(baker,
        main = "Baker Network Data",
        mar = c(1, 1, 3, 1),
        title.line = 2)

# SECTION 4.1.: Binary Blockmodeling --------------------------------------

# Data
bakerBinar <- baker
bakerBinar[bakerBinar > 0] <- 1

## Structural equivalence
resBinStr <- optRandomParC(
  M = bakerBinar,
  k = 3,
  rep = 1000,
  nCores = 0,
  blocks = c("nul", "com"),
  approach = "bin"
)

# The image matrix
IM(resBinStr)

# Visualization of the The image matrix
plot(
  resBinStr,
  main = "A Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)

# Matrix representation of the network of journals partitioned
# into 3 clusters using binary blockmodeling with structural
# equivalence and the corresponding block densities.
par(mfrow=c(1,2))
plotMat(
  M = bakerBinar,
  clu = clu(resBinStr),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(
  M = funByBlocks(resBinStr),
  main = "Block densities",
  mar = c(1, 1, 3, 1),
  title.line = 1
)

## Regular equivalence
resBinReg <- optRandomParC(
  M = bakerBinar,
  k = 2,
  rep = 1000,
  nCores = 0,
  blocks = c("nul", "com", "reg"),
  approach = "bin"
)

# Simple plot
plot( resBinReg)
# The network of journals partitioned into 2 clusters
# using binary blockmodeling with regular equivalence
# and the corresponding block densities.
par(mfrow=c(1,2))
plotMat(
  M = bakerBinar,
  clu = clu(resBinReg),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(
  M = funByBlocks(resBinReg),
  main = "Block densities",
  mar = c(1, 1, 3, 1),
  title.line = 1
)

# Adjusted Rand Index
crand2(clu1 = clu(resBinStr), clu2 = clu(resBinReg))

# SECTION 4.2.: Valued blockmodeling --------------------------------------

# Distribution of the number of citations among the journals.
hist(
  baker[baker > 0],
  breaks = 33,
  main = "",
  las = 1,
  xlab = "Number of citations"
)

# Summary statistics for non-zero values
summary(baker[baker > 0])

## Structural equivalence
resValStr <- optRandomParC(
  M = baker,
  k = 3,
  rep = 1000,
  preSpecM = 13,
  approach = "val",
  blocks = c("nul", "com"),
  nCores = 0
)

# simple plot
plot(resValStr)

# The network of journals partitioned into 3 clusters
# using valued blockmodeling (m=13) with structural
# equivalence and the corresponding block densities.
plotMat(
  M = baker,
  clu = clu(resValStr),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(
  M = funByBlocks(resValStr),
  main = "Block densities",
  mar = c(1, 1, 3, 1),
  title.line = 1
)

## Regular equivalence
resValReg <- optRandomParC(
  M = baker,
  k = 2,
  rep = 1000,
  preSpecM = 13,
  approach = "val",
  blocks = c("nul", "com", "reg"),
  nCores = 0,
  regFun = "max"
)

## Regular equivalence
resValRegK3 <- optRandomParC(
  M = baker,
  k = 3,
  rep = 1000,
  preSpecM = 13,
  approach = "val",
  blocks = c("nul", "com", "reg"),
  nCores = 0,
  regFun = "max"
)
plot(resValRegK3)
plot(resValRegK3,2) #second partition
canClu(partitions(resValRegK3)) #make all partitions canonical

# The network of journals partitioned into 2 clusters 
# using valued blockmodeling (m=13) with max regular 
# equivalence and the corresponding block densities.
plotMat(
  M = baker,
  clu = clu(resValReg),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(
  M = funByBlocks(resValReg),
  main = "Block densities",
  mar = c(1, 1, 3, 1),
  title.line = 1
)

# SECTION 4.3.: Homogenity blockmodeling  ---------------------------------

## Structural equivalence 
resHomSSStr <- optRandomParC(
  M = baker,
  k = 2,
  rep = 1000,
  approach = "hom",
  homFun = "ss",
  blocks = c("nul", "com"),
  nCores = 0
)

# The network of journals partitioned into 3 clusters
# using homogeneity blockmodeling (sum of squares) 
# structural equivalence and the corresponding block densities.
plotMat(
  M = baker,
  clu = clu(resHomSSStr),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(M = funByBlocks(resHomSSStr),
        main = "Block densities",
        mar = c(1, 1, 3, 1),
        title.line = 1)


## Regular equivalence
resHomSSReg <- optRandomParC(
  M = baker,
  k = 2,
  rep = 1000,
  approach = "hom",
  blocks = c("nul", "com", "reg"),
  regFun = "max",
  nCores = 0
)
plot(resHomSSReg)
IM(resHomSSReg)
# SECTION 4.4.: Pre-specified blockmodeling -------------------------------

# Image matrix
preImageReg <- rbind(c("com", "reg"),
                     c("reg", "nul"))

# Image matrix must be specified as an array
preImageRegCom <- array(NA, dim = c(2, 2, 2))
preImageRegCom[1, , ] <- rbind(c("com", "reg"),
                               c("reg", "nul"))
preImageRegCom[2, , ] <- rbind(c(NA, "com"),
                               c("com", "nul"))


# pre-specified blockmodeling
resValPre <- optRandomParC(
  M = baker,
  k = 2,
  rep = 1000,
  preSpecM = 13,
  approach = "val",
  blocks = preImageRegCom,
  nCores = 0
)

# The obtained image matrix
IM(resValPre)
plot(resValPre)

# The network of journals partitioned into 2 clusters
# using pre-specified valued blockmodeling and the corresponding block densities.
plotMat(
  M = baker,
  clu = clu(resValPre),
  main = "Baker Network Data",
  mar = c(1, 2, 3, 1),
  title.line = 2
)
plotMat(
  M = funByBlocks(resValPre),
  main = "Block densities",
  mar = c(1, 1, 3, 1),
  title.line = 1
)

