# ABOUT -------------------------------------------------------------------
# This is supplementary material.
# Article title: Scientific collaboration of researchers and organizations: A two-level blockmodeling approach
# Journal name: Scientometrics
# Authors: Marjan Cugmas (marjan.cugmas@fdv.uni-lj.si), Franc Mali (franc.mali@fdv.uni-lj.si) and Aleš Žiberna (ales.ziberna@fdv.uni-lj.si)
# Affiliation: Faculty of Social Sciences, University of Ljubljana, Kardeljeva ploščad 5, SI-1000 Ljubljana

# INSTALLING AND LOADING PACKAGES -----------------------------------------

install.packages("kmBlock", repos = "http://R-Forge.R-project.org")
library(kmBlock)
library(doRNG)
library(blockmodeling)
library(igraph)
library(randomcoloR)
library(car)

# SUPPORTING FUNCTIONS -----------------------------------------------------

source("supportingFunctions.R")
source("repairIgraphPlot.R")

# IMPORTING NETWORK DATA --------------------------------------------------

# The object network data is a list with the following elements:
# --- network The network is in a matrix form. The individuals and the organizations 
# are by lines and columns. The 1s indicate a link between two 
# individuals, two organizations or between individual and organization.
# --- clu A vector which defines if a given node is an individual 
# or an organization (1=individual level; 2=organizational level).
# --- fieldResearchers A vector: scientific field of an individual.
# 1 - Educational studies
# 2 - Economics
# 3 - Sociology
# 4 - Administrative and organizational sciences
# 5 - Criminology and social work
# 5 - Information science and librarianship
# 5 - Architecture and design
# 5 - Ethnic studies
# 6 - Political science
# 9 - Psychology
# 10 - Sport
# --- fieldOrganisatons A vector: scientific field of an organisation.
# 1 - Natural and biotechnical sciences
# 2 - Engineering sciences and technologies
# 3 - Medical sciences
# 4 - General social sciences
# 5 - Other social sciences
# 6 - Economics
# 7 - Humanities

# load the network data and the partition defining the level 
allF <- readRDS("networkData.RDS")
# extract the network
allFNet <- allF$network
# count the number of individuals and the number of organizations
allN <- unclass(table(allF$clu))

# CALCULATE THE WEIGHTS ---------------------------------------------------

# calculate the number of errors that would have been obtained if 
# each part would consists of just one block
EM <- blockmodeling::critFunC(allFNet,
                              clu = rep(1:2,times = allN),
                              approaches = "hom", 
                              blocks = "com")$EM[1,,]
# consider the inversively proportional values
wFix <- EM[1,1]/EM
# set the weight 0 to the empty part of the network
wFix[1,2] <- 0
# devide the weights for the two-mode part by two
wFix_1_2 <- wFix
wFix_1_2[2,1] <- wFix[2,1]/2
# apply these weights to all the network values/cells
cellWeights_1_2 <- expandMat(wFix_1_2, nn = allN)

# CREATING DIRECTORIES FOR RESULTS ----------------------------------------

dir.create("BMresults")
dir.create("BMfigs")
dir.create("GRAPHfig")

# APPLY BLOCKMODELING -----------------------------------------------------

# set the number of clusters for individual level 
kAut <- 2:25
# set the number of clusters for organizational level 
kOrg <- 2:8
# exand grid (apply blockmodeling for all combinations)
allCombinations <- expand.grid(kAut, kOrg)
# apply blockmodeling
# Applying blockmodeling is computationally very intensive and 
# it can take a considerable amount of time. 
# Change the above FALSE to TRUE to run the code.
if (FALSE) {
  for (i in 1:nrow(allCombinations)) {
    resAllMLk <- kmBlockORPC(M = allFNet, 
                             n = allN, 
                             k = unlist(allCombinations[i,]), 
                             rep = 1000, 
                             nCores = 0, 
                             weights = cellWeights_1_2)
    
    # save the results for later analysis
    saveRDS(resAllMLk, file = paste0("BMresults/BMres_", allCombinations[i, 1], "_", allCombinations[i, 2], ".RDS"))
  }
}

# VALUES OF THE CRITERION FUNCTION -------------------------------------

allCombinationsCF <- cbind(allCombinations, rep(NA, nrow(allCombinations)))
colnames(allCombinationsCF) <- c("k (aut)", "k (org)", "CF")
for (i in 1:nrow(allCombinations)) {
  BMres <- readRDS(paste0("BMresults/BMres_", allCombinations[i, 1], "_", allCombinations[i, 2], ".RDS"))
  allCombinationsCF[i,3] <- err(BMres)
}
gridPlotERR(data = allCombinationsCF)

# VISUALISATION - A NETWORK IN A MATRIX FORM ------------------------------

for (i in 1:nrow(allCombinations)) {
  BMres <- readRDS(paste0("BMresults/BMres_", allCombinations[i, 1], "_", allCombinations[i, 2], ".RDS"))
  cairo_pdf(file = paste0("BMfigs/BMfig_", allCombinations[i, 1], "_", allCombinations[i, 2], ".pdf"), width = 9, height = 8)
  plotMat(allFNet, 
          clu = clu(BMres),
          par.line.width = c(rep(0.1, allCombinations[i,1] - 1), 2, rep(0.1, allCombinations[i,2])),
          par.line.col = c(rep("blue", allCombinations[i,1] - 1), "#00000080", rep("red", allCombinations[i,2])))
  dev.off()
}

# VISUALISATION - A NETWORK IN A GRAPH FORM -------------------------------

for (i in 1:nrow(allCombinations)) {
  # load the blockmodeling result
  BMres <- readRDS(paste0("BMresults/BMres_", allCombinations[i,1], "_", allCombinations[i,2], ".RDS"))
  
  # set the number of clusters for individual level 
  kAut <- BMres$initial.param$k[1]
  # set the number of clusters for organizational level 
  kOrg <- BMres$initial.param$k[2]
  
  # determine vertex types
  vertexType <- rep(c(1, 0), BMres$initial.param$k)
  
  # calculate vertex sizes
  vertexSize <- table(clu(BMres))
  vertexSize[1:kAut] <- sqrt(vertexSize[1:kAut]) + 5
  nEmployByOrg <- rowSums(BMres$M[(BMres$initial.param$n[1] + 1):sum(BMres$initial.param$n), 1:BMres$initial.param$n[1]])
  belongingOrgClu <- clu(BMres)[(BMres$initial.param$n[1] + 1):sum(BMres$initial.param$n)]
  vertexSize[(kAut + 1):(kAut + kOrg)] <- sqrt(as.vector(by(nEmployByOrg, belongingOrgClu, sum))) + 5
  
  # image matrix - transofrm values for edge sizes
  IM <- getIM(BMres)
  resOmr <- sqrt(IM[vertexType == 1, vertexType == 1])
  IM[vertexType == 1, vertexType == 1] <- resOmr
  orgOmr <- sqrt(IM[vertexType == 0, vertexType == 0])
  IM[vertexType == 0, vertexType == 0] <- orgOmr
  freqResOrg <- fun.by.blocks(BMres$M, clu = clu(BMres), FUN = sum)[vertexType == 0,  vertexType == 1]
  relFreqResOrg <- freqResOrg/matrix(rep(table(clu(BMres))[vertexType == 1], kOrg), nrow = kOrg, byrow = TRUE)
  IM[vertexType == 0,  vertexType == 1] <- relFreqResOrg
  IM[upper.tri(IM)] = t(IM)[upper.tri(IM)]
  
  # define a network 
  network <- graph_from_adjacency_matrix(IM, weighted = TRUE, mode = "undirected", diag = TRUE)
  V(network)$vertexType <- vertexType
  V(network)$vertexSize <- vertexSize
  
  # define the order of drawing edges
  G2 <- network
  df <- as_data_frame(network)
  orderEdges <- order(df$weight)
  df2 <- df[orderEdges,]
  orderEdges <- as.numeric(unique(df2$from)) + 1
  G2 <- graph_from_data_frame(df2, directed = FALSE)
  
  # calculate values for drawing pies
  values <- by(c(allF[["fieldResearchers"]], allF[["fieldOrganisatons"]]), clu(BMres), table)
  myColor <- randomcoloR::distinctColorPalette(k = 10)
  
  # colors of the edges
  dfColors <- df2
  dfColors[,1] <- as.numeric(dfColors$from) + 1
  dfColors[,2] <- as.numeric(dfColors$to) + 1
  dfColors$colors <- rep(NA, nrow(dfColors))
  dfColors$bw <- sqrt(dfColors$weight)
  oo <- (dfColors$from %in% which(vertexType == 0)) & (dfColors$to %in% which(vertexType == 0))
  ll <- (dfColors$from %in% which(vertexType == 1)) & (dfColors$to %in% which(vertexType == 1))
  ol <- (dfColors$from %in% which(vertexType == 0)) & (dfColors$to %in% which(vertexType == 1))
  lo <- (dfColors$from %in% which(vertexType == 1)) & (dfColors$to %in% which(vertexType == 0))
  dfColors$colors[oo] <- rgb(red = 1, green = 1 - dfColors[oo, "bw"], blue = 1 - dfColors[oo, "bw"], maxColorValue = 1)
  dfColors$colors[ll] <- rgb(red = 1 - dfColors[ll, "bw"], green = 1 - dfColors[ll, "bw"], blue = 1, maxColorValue = 1)
  dfColors$colors[ol] <- rgb(red = 1 - dfColors[ol, "bw"], green = 1 - dfColors[ol, "bw"], blue = 1 - dfColors[ol, "bw"], maxColorValue = 1)
  dfColors$colors[lo] <- rgb(red = 1 - dfColors[lo, "bw"], green = 1 - dfColors[lo, "bw"], blue = 1 - dfColors[lo, "bw"], maxColorValue = 1)
  
  # plot the network
  cairo_pdf(file = paste0("GRAPHfig/fig_", paste(c(kAut, kOrg), collapse = "_"), ".pdf"), width = 9, height = 9)
  par(lwd = 2, mar = rep(0, 4), bg = "white")
  plot.igraph(G2,  
              loop.size = 0.3,
              vertex.size = vertexSize[orderEdges],
              vertex.shape = "pie", 
              vertex.pie = values[orderEdges],
              vertex.pie.color = list(myColor),
              pie.density = 5,
              pie.lty = 2,
              vertex.frame.color = rep(c("blue","red"), times = c(kAut, kOrg))[orderEdges], 
              vertex.label.dist = vertexSize[orderEdges]*0.1,
              vertex.label.color = "black",
              vertex.label.cex = 1.5,
              vertex.label = c(paste0("R", 1:kAut), paste0("O", 1:kOrg))[orderEdges],
              edge.width = E(G2)$weight*10,                     
              edge.arrow.size = 1,                          
              edge.arrow.width = 1,            
              edge.color = dfColors$colors,
              edge.loop.angle = -2,
              edge.lty = c("solid")
  )
  dev.off()
}