##### tex SAR network
# https://rdrr.io/cran/network/man/emon.html
# Drabek et al. (1981) provide seven case studies of emergent multi-organizational networks (EMONs) in the context of search and rescue (SAR) activities. Networks of interaction frequency are reported, along with several organizational attributes. 
# Network was recoded so that higher value means more frequent communication -> 1->4, 2->3, 3->2, 4->1 
# The tie values now mean:
# 0 - no communication
# 1 - “about once a day or less
# 2 - “every few hours”
# 3 - “about once an hour”
# 4 - “continuously”

library(blockmodeling)
library(sna)


# DATA --------------------------------------------------------------------

# Loading data
tex <- loadmatrix("texas.net") #File is in Pajek matrix format.
texCounty <- loadvector("texasCounty.clu")
texCounty <- factor(texCounty, labels = c("Other","Kendal","Kerr","Bandera"))

## Plotting the network and partition in matrix and graph form
plotMat(tex, clu = texCounty)
gplot1(tex/2, vertex.col = texCounty)

## Plotting the partitioned network and correspoding "image" (means by blocks)
par(mfrow=c(1,2))
plotMat(tex, 
        clu = texCounty, 
        mar=c(1,4,7,1),
        x.axis.val.pos = 1.02, 
        y.axis.val.pos = -0.02, 
        main = "Matrix by counties", title.line = 6,
        cex.val = 1)

tmpDensIM <- funByBlocks(tex, clu = texCounty, na.rm = TRUE)
plotMat(tmpDensIM,
        title.line = 2,
        main = "Density (image) matrix\ncounties",
        mar = c(1, 4, 7, 1),
        x.axis.val.pos = 1.03, 
        y.axis.val.pos = -0.03)
par(mfrow=c(1,1))


# HOMOGENEITY BLOCKMODELING -----------------------------------------------

## Homogeneity blockmodeling - structural equivalence
texStrSSk3 <- optRandomParC(M = tex, 
                            k = 3, 
                            rep = 200, 
                            blocks = c("com"), 
                            approach = "hom", 
                            homFun = "ss", 
                            nCores = 0)

texStrSSk3
plot(texStrSSk3, cex.val = 0.5)
plotMat(funByBlocks(texStrSSk3))

## Making a nice image matrix
strSSK3im <- funByBlocks(texStrSSk3, na.rm = TRUE, orderClu = TRUE)
rownames(strSSK3im) <- colnames(strSSK3im) <- c("Core","Kerr","Pheriphery")
plotMat(strSSK3im,
        main = "Image matrix", 
        mar = c(1, 4, 7, 1), 
        x.axis.val.pos = 1.03, 
        y.axis.val.pos = -0.03, 
        title.line = 6)

## Selecting number of clusters based on the value of a criterion function - using a for loop
texStrSSList<-list()
for(k in 1:5){
  texStrSSList[[as.character(k)]]<-optRandomParC(M = tex,
                                                 k = k,
                                                 rep = ifelse(k==1,1,200),
                                                 blocks = c("com"),
                                                 approach = "hom",
                                                 homFun = "ss",
                                                 nCores = ifelse(k==1,1,0))
}
plot(sapply(texStrSSList,err)~I(1:5),type="o", xlab="k",ylab="CF")

## Selecting number of clusters based on the Relative fit 
# https://link.springer.com/article/10.1007/s10260-021-00595-1 
# The code below can be time consuming
if(file.exists("texStrSSRF.Rdata")){
    load("texStrSSRF.Rdata")  
  }else{
	set.seed(1)
	texStrSSRF <- NULL
	for(k in 1:8){
	  tmp <- optRandomParC(M = tex, 
								  k = k, 
								  rep = 200, 
								  blocks = c("com"), 
								  approach = "hom", 
								  homFun = "ss", 
								  nCores = 0)
	  
	  # m stands for the number of randomized networks and
	  # should be increased to, e.g., 30
	  texStrSSRF[k] <- RF(tmp, m = 5)$RF
	}
	save(texStrSSRF,file = "texStrSSRF.Rdata")
}
plot(y = texStrSSRF, x = 1:length(texStrSSRF), type="o", xlab = "k", ylab = "RF")

# BINARY BLOCKMODELING ----------------------------------------------------

## structural equivalence
texStrBink3 <- optRandomParC(M = tex>0, 
                             k = 3, 
                             rep = 200, 
                             blocks = c("nul", "com"), 
                             approach = "bin", 
                             nCores = 0)
texStrBink3
plot(texStrBink3)
plot(texStrBink3, orderClu=TRUE)

# image matrix
plotMat(funByBlocks(texStrBink3))
plotMat(funByBlocks(texStrBink3, orderClu = TRUE))

gplot1(funByBlocks(texStrBink3))
gplot1(tex>0,vertex.col = clu(texStrBink3))

## comparison with homogeneity partition
table(clu(texStrSSk3),clu(texStrBink3))

## regular equivalence
# 3 clusters
set.seed(5)
texRegBink3 <- optRandomParC(M = tex>0,
                             k = 3,
                             rep = 200,
                             blocks = c("nul","com", "reg"),
                             approach = "bin",
                             nCores = 0)
texRegBink3
plot(texRegBink3)
plot(texRegBink3, which = 2)
table(clu(texRegBink3), clu(texRegBink3, 2))

# 2 clusters
texRegBink2<-optRandomParC(M = tex>0, 
                           k = 2, 
                           rep = 200, 
                           blocks=c("nul", "com", "reg"), 
                           approach = "bin", 
                           nCores = 0)
plot(texRegBink2)
texRegBink2

# VALUED BLOCKMODELING ----------------------------------------------------

## Structural equivalence, m = 4
table(tex)
summary(tex[tex>0])
set.seed(2021)
texStrValm4k3<-optRandomParC(M = tex, 
                             preSpecM = 4, 
                             k = 3, 
                             rep = 200, 
                             blocks=c("nul", "com"), 
                             approach = "val", 
                             nCores = 0)
plot(texStrValm4k3)
texStrValm4k3			 
err(texStrValm4k3)
funByBlocks(texStrValm4k3)
table(clu(texStrValm4k3),texCounty)								   

texStrValm4k3im<-funByBlocks(texStrValm4k3, orderClu = TRUE)
rownames(texStrValm4k3im) <- colnames(texStrValm4k3im)<-c("Core", "Kerr", "Pheriphery")
plotMat(texStrValm4k3im)

## Structural equivalence, m = 2
table(tex)
summary(tex[tex>0])
set.seed(2021)
texStrValm2k3 <- optRandomParC(M = tex, 
                               preSpecM = 2, 
                               k = 3, 
                               rep = 200, 
                               blocks = c("nul", "com"), 
                               approach = "val", 
                               nCores = 0)
texStrValm2k3
plot(texStrValm2k3)
funByBlocks(texStrValm2k3)

## Regular equivalence, m = 4 ("does not work")
set.seed(2021)
texRegValm4k3 <- optRandomParC(M = tex,
                               preSpecM = 4,
                               k = 3,
                               rep = 200,
                               blocks = c("nul", "com", "reg"),
                               approach = "val", 
                               nCores = 0)
plot(texRegValm4k3)

# PRE-SPECIFIED BLOCKMODELING ---------------------------------------------

## Cohesive groups
IMcoh2 <- ifelse(diag(2)==1, "com","nul")
IMcoh3 <- ifelse(diag(3)==1, "com","nul")
IMcoh4 <- ifelse(diag(4)==1, "com","nul")
IMcoh4

## errors for one and "maximal" (same as number of units) cohesive groups
c(sum(4-tex)-4*nrow(tex), sum(tex))
set.seed(2021)
texStrValm4coh2<-optRandomParC(M = tex,
                               preSpecM = 4,
                               k = 2,
                               rep = 200,
                               blocks = IMcoh2,
                               approach = "val",
                               nCores = 0)
table(clu(texStrValm4coh2), clu(texStrValm4coh2,which = 2))

set.seed(2021)
texStrValm4coh2 <- optRandomParC(M = tex,
                                 preSpecM = 4,
                                 k = 2,
                                 rep = 200,
                                 blocks = IMcoh2,
                                 approach = "val",
                                 nCores = 0, 
                                 switch.names = TRUE)
plot(texStrValm4coh2)
plot(texStrValm4coh2, which = 2)
table(clu(texStrValm4coh2), clu(texStrValm4coh2, which = 2))

texStrValm4coh3 <-optRandomParC(M = tex,
                                preSpecM = 4,
                                k = 3,
                                rep = 200,
                                blocks = IMcoh3,
                                approach = "val",
                                nCores = 0,
                                switch.names = TRUE)
plot(texStrValm4coh3)

texStrValm4coh4 <- optRandomParC(M = tex,
                                  preSpecM = 4,
                                  k = 4,
                                  rep = 200,
                                  blocks = IMcoh4,
                                  approach = "val",
                                  nCores = 0,
                                  switch.names = TRUE)
plot(texStrValm4coh4)

## making a vector or errors
errCoh<-c(sum(4-tex)-4*nrow(tex), err(texStrValm4coh2), err(texStrValm4coh3), err(texStrValm4coh4))
plot(errCoh,type="o")
abline(h=sum(tex))

if(FALSE){
  if(file.exists("texStrValm4cohList.Rdata")){
    load("texStrValm4cohList.Rdata")  
  }else{
    texStrValm4cohList<-list()
    for(k in 2:15){
      tmpIM<-ifelse(diag(k)==1, "com","nul")
      texStrValm4cohList[[as.character(k)]]<-optRandomParC(M=tex,
                                                           preSpecM=4,
                                                           k=k,
                                                           rep=100,
                                                           blocks=tmpIM,
                                                           approach="val",
                                                           nCores=0,
														   switch.names = TRUE)
    }
    save(texStrValm4cohList,file = "texStrValm4cohList.Rdata")
  }
  
  errCoh<-c("1"=sum(4-tex), sapply(texStrValm4cohList,err))
  plot(errCoh,type="o")
  abline(h=sum(tex))
  plot(errCoh,type="o",ylim=c(min(errCoh),650))
  abline(h=sum(tex))
}


## Core-perhiphery - structural equivalence
# model:
#      1        2
# 1   com    nul,com
# 2 nul,com    nul
IMcpStr <- array(NA, dim=c(2,2,2))
IMcpStr[1,,] <- "nul"
IMcpStr[1,1,1] <- "com"
IMcpStr[2,2,1] <- "com"
IMcpStr[2,1,2] <- "com"
IMcpStr<-array(NA,dim=c(2,2,2))
IMcpStr[,1,1]<-c("com", NA)
IMcpStr[,2,2]<-c("nul", NA)
IMcpStr[,1,2]<-c("nul", "com")
IMcpStr[,2,1]<-c("nul", "com")

# for better idea what we have done
printBlocks(IMcpStr)

texStrValm4CPstr <- optRandomParC(M = tex,
                                  preSpecM = 4,
                                  k = 2,
                                  rep = 200,
                                  blocks = IMcpStr,
                                  approach = "val",
                                  nCores = 0)
plot(texStrValm4CPstr)
texStrValm4CPstr


## Core-perhiphery - regular equivalence
# model:
#        1              2
# 1     com        nul,com,reg
# 2 nul,com,reg        nul
IMcpReg<-array(NA,dim=c(3,2,2))
IMcpReg[1,,]<-"nul"
IMcpReg[1,1,1]<-"com"
IMcpReg[2,2,1]<-"com"
IMcpReg[2,1,2]<-"com"
IMcpReg[3,2,1]<-"reg"
IMcpReg[3,1,2]<-"reg"

# for better idea what we have done
printBlocks(IMcpReg)

texStrValm4CPreg <- optRandomParC(M = tex,
                                  preSpecM = 4,
                                  k = 2,
                                  rep = 200,
                                  blocks = IMcpReg,
                                  approach = "val",
                                  nCores = 0)
plot(texStrValm4CPreg)
err(texStrValm4CPreg)
plot(texStrValm4CPreg, use.IM=FALSE)

# Relative Fit can be used to compare different blockmodels
# The code below can be time consuming
if(file.exists("RFvec.Rdata")){
	load("RFvec.Rdata")
}else {
  set.seed(1)
	RFvec<-c(
		texStrValm4coh3 = RF(texStrValm4coh3)$RF,
		texStrValm4coh4 = RF(texStrValm4coh4)$RF,
		texStrValm4CPreg = RF(texStrValm4CPreg)$RF
	)
	save(RFvec, file = "RFvec.Rdata")
}
barplot(RFvec)
##### Example of a more complex model
## Make a model for this situation
# - groups should be cohesive, except one, which does not need to be (but can be)
# - cohesive groups should be complete
# - ties between groups can be anything (nul, com, reg)

makeMyIM<-function(k){
  myIM<-array(NA,dim=c(3,k,k))
  myIM[1,,]<-"nul"
  myIM[cbind(1,1:(k-1),1:(k-1))]<-"com"
  myIM[2,,]<-"com"
  myIM[cbind(2,1:(k-1),1:(k-1))]<-NA
  myIM[3,,]<-"reg"
  myIM[cbind(3,1:k,1:k)]<-NA
  return(myIM)
}
printBlocks(makeMyIM(4))

k<-3; texStrValm4myIMk3<-optRandomParC(M = tex, preSpecM=4,k=k,rep=200,blocks=makeMyIM(3),approach="val",nCores=0)
plot(texStrValm4myIMk3)
printBlocks(IM(texStrValm4myIMk3))
plot(texStrValm4myIMk3,which = 2)
printBlocks(IM(texStrValm4myIMk3,which = 2))



