# Social network analysis package is needed to draw the network
require(sna)

GroupBasedRandDigraph <- function(Mat,Part){
  S <- dim(Mat)[1]
  # check whether the partition starts in 0 or 1
  if (min(Part)==0) Part <- Part+1
  # Find a permutation that order the nodes according to groups
  OrderPart<-sort(Part,index.return=T)$ix
  # Order the matrix
  OrderMat <- Mat[OrderPart,OrderPart]
  # Order the partition
  NewPart <- sort(Part)
  # K is the number of groups
  K <- max(NewPart)
  # Size sumbmatrices ij = size group i x size group j
  Sizeij <- outer(table(NewPart),table(NewPart))
  # Count Lij = number of links from nodes in i to nodes in j
  Agg <- matrix(0,S,S)
  for (i in 1:S) Agg[i,NewPart[i]] <- 1
  Lij <- t(Agg) %*% OrderMat %*% (Agg)
  Lij <- Lij[1:K,1:K]
  rownames(Lij)<-rownames(Sizeij,prefix="")
  colnames(Lij)<-rownames(Sizeij,prefix="")
  rownames(Sizeij) <- rownames(Lij,prefix="")
  colnames(Sizeij) <- colnames(Lij,prefix="")
  # Probabilities are obtained dividing Lij/Sizeij
  Pij <- Lij/Sizeij
  # Compute log-likelihoods
  LogLMat <- Lij*log(Pij)+(Sizeij-Lij)*log(1-Pij)
  # Remove NaN (that appear when pij is 0 or 1
  LogLMat[is.nan(LogLMat)] <- 0
  # Print output
  print("L_ij matrix")
  print(Lij)
  print("Size_ij matrix")
  print(Sizeij)
  print("P_ij matrix")
  print(Pij)
  print("LogLikelihood")  
  LogL <- sum(LogLMat)
  print(LogL)
  print("AIC")
  AIC <- 2*S+2*K*K-2*LogL
  print(AIC)
  OriginalNames <- 1:S
  OrderedNames <- OriginalNames[OrderPart]
  rownames(OrderMat) <- c(OrderedNames)
  colnames(OrderMat) <- rownames(OrderMat)
  # DRAW NETWORK
  Colo <- rainbow(K+2)
  layout(matrix(1:2,1,2))
  gplot(OrderMat,vertex.col=Colo[NewPart+1],displaylabels=T,label=OrderedNames)
  ColorMat <- matrix(0,S+1,S+1)
  ColorMat[1,2:(S+1)] <- NewPart+1
  ColorMat[2:(S+1),1] <- NewPart+1
  ColorMat[2:(S+1),2:(S+1)] <- OrderMat
  rownames(ColorMat)[2:(S+1)] <- rownames(OrderMat)
  colnames(ColorMat)[2:(S+1)] <- rownames(OrderMat)
  rownames(ColorMat)[1] <- "k"
  colnames(ColorMat)[1] <- "k"
  print(ColorMat)
  plot.mygroupmatrix(ColorMat,GroupColors=Colo)
  return(AIC)
}

plot.mygroupmatrix <-function (x, labels = NULL, drawlab = TRUE, diaglab = TRUE, drawlines = TRUE, xlab = NULL, ylab = NULL, cex.lab = 1, GroupColors,...) 
{
    if (class(x) == "network") 
        x <- as.sociomatrix.sna(x)
    n <- dim(x)[1]
    o <- dim(x)[2]
    if (is.null(labels)) 
        labels <- list(NULL, NULL)
    if (is.null(labels[[1]])) {
        if (is.null(rownames(x))) 
            labels[[1]] <- 1:dim(x)[1]
        else labels[[1]] <- rownames(x)
    }
    if (is.null(labels[[2]])) {
        if (is.null(colnames(x))) 
            labels[[2]] <- 1:dim(x)[2]
        else labels[[2]] <- colnames(x)
    }
    #d <- 1 - (x - min(x, na.rm = TRUE))/(max(x, na.rm = TRUE) - 
    #    min(x, na.rm = TRUE))
    d <- x
    if (is.null(xlab)) 
        xlab <- ""
    if (is.null(ylab)) 
        ylab <- ""
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n", 
        axes = FALSE, xlab = xlab, ylab = ylab, ...)
    for (i in 1:n){
      for (j in 1:o) {
        if (d[i,j]==0) rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = "white", xpd = TRUE, border = drawlines)
        if (d[i,j]==1)rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = "black", xpd = TRUE, border = drawlines)
        if (d[i,j]>1) rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = GroupColors[d[i,j]], xpd = TRUE, border = drawlines)
      }
    }
    rect(0.5, 0.5, o + 0.5, n + 0.5, col = NA, xpd = TRUE)
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]][1:(n)], cex = cex.lab)
        text(1:o, rep(0, o), labels[[2]][1:(n)], cex = cex.lab)
    }
    if ((n == o) & (drawlab) & (diaglab)) 
        if (all(labels[[1]] == labels[[2]])) 
            text(1:o, 1:n, labels[[1]][1:(n)], cex = cex.lab)
}