########################################
# John Rigsby
# RigsbyJT@nswc.navy.mil
# 
# X is an affiliation network matrix of size m by n where
# m is the number of rows and n is the number of columns
#
# Orx <- com.row(X) makes Orx the row overlap matrix of matrix X
# Ocx <- com.col(X) makes Ocx the column overlap matrix of matrix X
# 
# ROrx <- com.row.rel(X) makes ROrx the relative row overlap matrix 
#         of matrix X.  It is the values of the overlap matrix where
#         the values of each row is divided by that rows diagonal value
# ROcx <- com.col.rel(X) makes ROrx the relative col overlap matrix 
#         of matrix X.  It is the values of the overlap matrix where
#         the values of each row is divided by that rows diagonal value
#
# IROrx <- com.row.rel.imp(X) makes IROrx the improved relative row overlap 
#         matrix of matrix X.  It is the values of the overlap matrix where
#         the values of each row is divided by that rows(i) diagonal value and 
#         other row (j) diagonal value summed.
# IROcx <- com.col.rel.imp(X) makes IROrx the improved relative col overlap 
#         matrix of matrix X.  It is the values of the overlap matrix where
#         the values of each row is divided by that rows(i) diagonal value and 
#         other row (j) diagonal value summed.
#
com.row.rel.imp <- function(X)
{
   # X is m by n with m actors and n events
   m <- nrow(X)
   n <- ncol(X)
   
   Orow <- as.matrix(X) %*% as.matrix(t(X))
   holder <- Orow
   # Orow is the co-membership matrix
   for(i in 1:m)
      for(j in 1:m)
         Orow[i,j] <- 2*Orow[i,j]/(holder[i,i]+holder[j,j])
   # return this
   Orow
}

com.row.rel <- function(X)
{
   # X is m by n with m actors and n events
   m <- nrow(X)
   n <- ncol(X)
   
   Orow <- as.matrix(X) %*% as.matrix(t(X))
   holder <- Orow
   # Orow is the co-membership matrix
   for(i in 1:m)
      for(j in 1:m)
         Orow[i,j] <- Orow[i,j]/holder[i,i]
   # return this
   Orow
}
com.row    <- function(X)
{
   # X is m by n with m actors and n events
   
   # Orow is the co-membership matrix
   Orow <- as.matrix(X) %*% as.matrix(t(X))
}
com.col.rel.imp <- function(X)
{
   # X is m by n with m actors and n events
   m <- nrow(X)
   n <- ncol(X)
   
   Ocol <- as.matrix(t(X)) %*% as.matrix(X)
   holder <- Ocol
   # Ocol is the co-membership matrix
   for(i in 1:n)
      for(j in 1:n)
         Ocol[i,j] <- 2*Ocol[i,j]/(holder[i,i]+holder[j,j])
   # return this
   Ocol
}
com.col.rel    <- function(X)
{
   # X is m by n with m actors and n events
   m <- nrow(X)
   n <- ncol(X)
   
   Ocol <- as.matrix(t(X)) %*% as.matrix(X)
   holder <- Ocol
   # Ocol is the co-membership matrix
   for(i in 1:n)
      for(j in 1:n)
         Ocol[i,j] <- Ocol[i,j]/holder[i,i]
   # return this
   Ocol
}
com.col <- function(X)
{
   # X is m by n with m actors and n events
   
   # Ocol is the co-membership matrix
   Ocol <- as.matrix(t(X)) %*% as.matrix(X)
}

onezero <- function(X)
{
    Y <- X
    for(i in 1:(nrow(X)))
       for(j in 1:(ncol(X))) 
          if(X[i,j]>0) 
              Y[i,j]=1
    # RETURN Y
    Y
}
           
zeroone <- function(X)
{
    Y <- X
    for(i in 1:(nrow(X)))
       for(j in 1:(ncol(X))) 
          if(X[i,j]!=1) 
              Y[i,j]=0
    # RETURN Y
    Y
}
           
hypd <- function(X,N)
{
   # N is the total number of events that the actors
   # in the comembership matrix could have gone to

   if(nrow(X) != ncol(X)) stop("Comembership matrix must be square")
   # X is Y by Y comembership matrix
   # Z will be the return matrix
   Z <- X 
   Y <- nrow(X)
   
   for(i in 1:Y){
      for(j in 1:Y) {
         k <- X[i,i] # total number of events actor A went to
         m <- X[j,j] # total number of events actor B went to
         n <- N - m # total number of events B did not goto
         x <- X[i,j] # the number of events actors A & B share
         mean <- (k*m/N) # expected mean

         if( mean == x ) {
            Z[i,j]=1
         }
         if( mean > x ) {
            # X[i,j] is the number of events actors A & B have in common
            lowx <- 0:x # lower tail of dist            
            if((2*mean-x)>0) {
              upx <- (mean+(mean-x)):N # upper tail of dist
              Z[i,j] = sum(dhyper(lowx,m,n,k)) +sum(dhyper(upx,m,n,k))
            } else {Z[i,j] = sum(dhyper(lowx,m,n,k))}
         } 
         if(mean < x ) {
            upx <- x:N # upper tail of dist
            if((2*mean-x)>0) {
              lowx <- 0:(mean - (x-mean)) # lower tail of dist
              Z[i,j] = sum(dhyper(lowx,m,n,k)) +sum(dhyper(upx,m,n,k))
            } else {Z[i,j] = sum(dhyper(upx,m,n,k))}
         }
      }
   }

   # Return Z
   Z
}

rhypd <- function(X,N)
{
   # N is the total number of events that the actors
   # in the comembership matrix could have gone to

   if(nrow(X) != ncol(X)) stop("Comembership matrix must be square")
   # X is Y by Y comembership matrix
   # Z will be the return matrix
   Z <- X 
   Y <- nrow(X)
   
   for(i in 1:Y){
      for(j in 1:Y) {
         k <- X[i,i] # total number of events actor A went to
         m <- X[j,j] # total number of events actor B went to
         n <- N - m # total number of events B did not goto
         x <- X[i,j] # the number of events actors A & B share
         mean <- round(k*m/N) # expected mean

         if( mean == x ) {
            Z[i,j]=1
         }
         if( mean > x ) {
            # X[i,j] is the number of events actors A & B have in common
            lowx <- 0:x # lower tail of dist
            upx <- (mean+(mean-x)):N # upper tail of dist
            Z[i,j] = sum(dhyper(lowx,m,n,k)) + sum(dhyper(upx,m,n,k))
         } 
         if(mean < x ) {
            lowx <- 0:(mean - (x-mean)) # lower tail of dist
            upx <- x:N # upper tail of dist
            Z[i,j] = sum(dhyper(lowx,m,n,k)) + sum(dhyper(upx,m,n,k))
         }
      }
   }

   # Return Z
   Z
}

exphypd <- function(X,N)
{
   # N is the total number of events that the actors
   # in the comembership matrix could have gone to

   if(nrow(X) != ncol(X)) stop("Comembership matrix must be square")
   # X is Y by Y comembership matrix
   # Z will be the return matrix
   Z <- X 
   Y <- nrow(X)
   
   for(i in 1:Y){
      for(j in 1:Y) {
         k <- X[i,i] # total number of events actor A went to
         m <- X[j,j] # total number of events actor B went to
         mean <- (k*m/N) # expected mean
         Z[i,j] = mean
      }
   }

   # Return Z
   Z
}

rexphypd <- function(X,N)
{
   # N is the total number of events that the actors
   # in the comembership matrix could have gone to

   if(nrow(X) != ncol(X)) stop("Comembership matrix must be square")
   # X is Y by Y comembership matrix
   # Z will be the return matrix
   Z <- X 
   Y <- nrow(X)
   
   for(i in 1:Y){
      for(j in 1:Y) {
         k <- X[i,i] # total number of events actor A went to
         m <- X[j,j] # total number of events actor B went to
         mean <- (k*m/N) # expected mean
         Z[i,j] = round(mean)
      }
   }

   # Return Z
   Z
}

ldtest  <- function(X)
{
    dcount <- 0
    dtot <- 0
    ndcount <- 0
    ndtot <- 0
    Y <- X
    max <- 0 

    for(i in 1:(nrow(X)))
       for(j in 1:(ncol(X))) 
       if (is.na(X[i,j])) { } else
       {
          if(max < X[i,j]) {max <- X[i,j] }
       }
    if(max<1) {max <- 1}

    for(i in 1:(nrow(X)))
       for(j in 1:(ncol(X))) 
       {
          if(is.na( X[i,j])) {X[i,j] <- max }
          X[i,j] <- X[i,j]/max
       }

    for(i in 1:(nrow(X)))
       for(j in 1:(ncol(X))) 
       {
          if(is.na(X[i,j])) { X[i,j] <- 0 }
          if(i == j) {
              dcount <- dcount + 1 
              dtot <- dtot + X[i,j]
          } else {
              ndcount <- ndcount + 1 
              ndtot <- ndtot + X[i,j]
          } 
       }
    dmean <-  dtot/dcount          
    ndmean <-  ndtot/ndcount          


    returnvalue <- dmean/ndmean

    returnvalue
}

myplotblockmodel <- function (x, ...)
{
    oldpar <- par()
    n <- dim(x$blocked.data)[2]
    m <- stackcount(x$blocked.data)
    if (!is.null(x$plabels))
        plab <- x$plabels
    else plab <- (1:n)[x$order.vector]
    if (!is.null(x$glabels))
        glab <- x$glabels
    else glab <- 1:m
    print(glab)
    par(mfrow = c(floor(sqrt(m)), ceiling(m/floor(sqrt(m)))))
    if (m > 1)
        for (i in 1:m) {
            plot.sociomatrix(x$blocked.data[i, , ], labels = list(plab,
                plab), main = paste("Relation - ", glab[i]))
            for (j in 2:n) if (x$block.membership[j] != x$block.membership[j - 1])
                abline(v = j - 0.5, h = j - 0.5, lty = 1,col="red",lwd=3)
        }
    else {
        plot.sociomatrix(x$blocked.data, labels = list(plab, plab),
            main = paste("Relation - ", glab[1]))
        for (j in 2:n) if (x$block.membership[j] != x$block.membership[j - 1])
            abline(v = j - 0.5, h = j - 0.5, lty = 1,col="red",lwd=3)
    }
    par(oldpar)
}

           
mymatrixplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if(( 0.025 < d[i,j]) && ( d[i,j] < 0.05)) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(d[i, j],0,1), xpd = TRUE) 
      } else if ( 0.025 > d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else if ( 0.975 < d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(0,d[i, j],0), xpd = TRUE) 
      } else if ( ( 0.975 >  d[i,j]) && (0.95 < d[i,j]) ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}

mymatrixlowplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if(( 0.025 < d[i,j]) && ( d[i,j] < 0.05)) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(d[i, j],0,1), xpd = TRUE) 
      } else if ( 0.025 > d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}

mymatrixhighplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if ( 0.975 < d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(0,d[i, j],0), xpd = TRUE) 
      } else if ( ( 0.975 >  d[i,j]) && (0.95 < d[i,j]) ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}


mymatrixlowbfplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if(( 0.025/420 < d[i,j]) && ( d[i,j] < 0.05/420)) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(d[i, j],0,1), xpd = TRUE) 
      } else if ( 0.025/420 > d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}

mymatrixhighbfplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if ( 1- 0.025/420 < d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(0,d[i, j],0), xpd = TRUE) 
      } else if ( ( 1- 0.05/420 <  d[i,j]) && (1 - 0.025/420 > d[i,j]) ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}

mymatrixbfplot <- function (x, labels = list(seq(1:dim(x)[1]), seq(1:dim(x)[2])),
    drawlab = TRUE, diaglab = TRUE, ...)
{
    n <- dim(x)[1]
    o <- dim(x)[2]
    d <- 1 - (x - min(x))/(max(x) - min(x))
    plot(1, 1, xlim = c(0, o + 1), ylim = c(n + 1, 0), type = "n",
        axes = FALSE, ...)
    for (i in 1:n) for (j in 1:o) 
      if(( 0.025/420 < d[i,j]) && ( d[i,j] < 0.05/420)) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(d[i, j],0,1), xpd = TRUE) 
      } else if ( 0.025/420 > d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else if ( 1- 0.025/420 < d[i,j] ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(0,d[i, j],0), xpd = TRUE) 
      } else if ( ( 1- 0.05/420 <  d[i,j]) && (1 - 0.025/420 > d[i,j]) ) {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = rgb(1,d[i, j],0), xpd = TRUE) 
      } else {
        rect(j - 0.5, i + 0.5, j + 0.5, i - 0.5, col = gray(d[i, j]), xpd = TRUE) 
      }
    if (drawlab) {
        text(rep(0, n), 1:n, labels[[1]])
        text(1:o, rep(0, o), labels[[2]])
    }
    if ((n == o) & (drawlab) & (diaglab))
        if (labels[[1]] == labels[[2]])
            text(1:o, 1:n, labels[[1]])
}

allegiance<-function(dat,ec,k=NULL,h=NULL,block.content="density",plabels=ec$plabels,glabels=ec$glabels,rlabels=NULL,mode="digraph"){
   #require(mva)
   #First, extract the blocks
   b<-cutree(ec$cluster,k,h)
   #Prepare the data
   n<-dim(dat)[2]
   if(length(dim(dat))>2)
      d<-dat
   else{
      d<-array(dim=c(1,n,n))
      d[1,,]<-dat
   }
#   if(!diag)
#      d<-diag.remove(d)
   #Now, construct a model

   rm <- dim(d)[1]
#   allval <- array(data = 0,dim = c(3,n))
   allval <- array(data = 0,dim = c(rm,3,n))
   blockhelp <- array(data=0,dim=c(k))
   blockhurt <- array(data=0,dim=c(k))
   blocktotal <- array(data=0,dim=c(k))
   blockratio <- array(data=0,dim=c(k))
   bnum <- array(data=0,dim=c(k))
   for(i in 1:rm){
      for (ac in 1:n) {
         bnum[b[ac]] = bnum[b[ac]]+1
         for(j in 1:n) {
            if(b[j]==b[ac]){
                allval[i,1,ac]=allval[i,1,ac]+d[i,ac,j]
#                allval[1,ac]=allval[1,ac]+d[i,ac,j]
            } else {
                allval[i,2,ac]=allval[i,2,ac]+d[i,ac,j]
#                allval[2,ac]=allval[2,ac]+d[i,ac,j]
            }
         }
      allval[i,2,ac] = allval[i,2,ac]
      }
   }
   for(x in 1:rm){
      for(blocknum in 1:k){
         for(y in 1:n){
            for(z in 1:n){
                if(b[y]==blocknum && b[z]==blocknum) {
                    blockhelp[blocknum]=blockhelp[blocknum]+d[x,y,z]
                } else if (b[y]!=blocknum && b[z]==blocknum) {
                    blockhurt[blocknum]=blockhurt[blocknum]+d[x,y,z]
                } else if (b[y]==blocknum && b[z]!=blocknum) {
                    blockhurt[blocknum]=blockhurt[blocknum]+d[x,y,z]
                }
            }
         }
      }
   }


   totalhelp = 0
   totalhurt = 0
   totalall = 0
   totalratio = 0
   for (ac in 1:n) {
#       allval[1,1,ac]=(allval[1,1,ac])/(bnum[b[ac]])
#       if(k==1) {
#              allval[i,2,ac]=0
#       } else {
#              allval[1,2,ac]=(allval[1,2,ac])/(n-bnum[b[ac]])
#       }
       allval[1,3,ac]=allval[1,1,ac]-allval[1,2,ac]
       totalhelp = totalhelp + allval[1,1,ac]
       totalhurt = totalhurt + allval[1,2,ac]
       totalall = totalall + allval[1,3,ac]
   }
   totalratio = totalhurt/totalhelp
   for (j in 1:k) {
         blockhelp[j] = blockhelp[j] /(bnum[j]*bnum[j])
         if(k!=1) {
             blockhurt[j] = blockhurt[j] /(2*bnum[j]*(n-bnum[j]))
             blocktotal[j] = blockhelp[j]-blockhurt[j]
             blockratio[j] = blockhurt[j]/blockhelp[j]
         } else {
             blockhurt[j] = 0
             blocktotal[j] = blockhelp[j]
             blockratio[j] = 0
         }
   }
         

   o<-list()
   o$block.membership<-b[ec$cluster$order]
   o$order.vector<-ec$cluster$order
   if(length(dim(allval))>2){
      o$blocked.data<-allval[1,,ec$cluster$order]
   }else{
      o$blocked.data<-allval[1,,ec$cluster$order]
   }
   if(dim(allval)[1]==1){
      o$allegiance.model<-allval[1,,]
      colnames(o$allegiance.model) = ec$glabels[ec$cluster$order]
      rownames(o$allegiance.model) = c("help", "hurt","allegiance")
   }else{
      o$allegiance.model<-allval
      dimnames(o$block.model)<-list(glabels,ec$glabels[ec$cluster$order],c("help", "hurt","allegiance"))
   }
   o$plabels<-plabels[ec$cluster$order]
 #  o$glabels<-glabels
   o$rlabels<-rlabels
   o$cluster.method<-ec$cluster.method
   o$equiv.fun<-ec$equiv.fun
   o$equiv.metric<-ec$metric
   o$total.allegiance.help<- totalhelp
   o$total.allegiance.hurt<- totalhurt
   o$total.allegiance.total<- totalall
   o$total.allegiance.ratio<- totalratio
   o$block.help <- blockhelp
   o$block.hurt <- blockhurt
   o$block.total <- blocktotal
   o$block.ratio <- blockratio
   o$block.num <- bnum



   class(o)<-"allegiance"
   o
}



myplotallegiance <- function (x, ...)
{
    oldpar <- par()
    n <- dim(x$blocked.data)[1]
    m <- stackcount(x$blocked.data)
    if (!is.null(x$plabels)){
        plab <- x$plabels
    }else {plab <- (1:n)[x$order.vector]}
    if (!is.null(x$glabels))
        glab <- x$glabels
    else glab <- 1:m
    print(glab)
    par(mfrow = c(floor(sqrt(m)), ceiling(m/floor(sqrt(m)))))
    if (m > 1)
        for (i in 1:m) {
            plot.sociomatrix(x$blocked.data[i, , ], labels = list(plab,
                plab), main = paste("Relation - ", glab[i]))
            for (j in 2:n) if (x$block.membership[j] != x$block.membership[j - 1])
                abline(v = j - 0.5, lty = 1,col="red",lwd=3)
        }
    else {
        plot.sociomatrix(x$blocked.data, labels = list(plab, plab),
            main = paste("Relation - ", glab[1]))
        for (j in 2:n) if (x$block.membership[j] != x$block.membership[j - 1])
            abline(v = j - 0.5, lty = 1,col="red",lwd=3)
    }
    par(oldpar)
}

           

newallegiance<-function(dat,ec,k=NULL,h=NULL,block.content="density",plabels=ec$plabels,glabels=ec$glabels,rlabels=NULL,mode="digraph"){
   #require(mva)
   #First, extract the blocks
   b<-cutree(ec$cluster,k,h)
   #Prepare the data
   n<-dim(dat)[2]
   if(length(dim(dat))>2){
      d<-dat
   }else{
      d<-array(dim=c(1,n,n))
      d[1,,]<-dat
   }
#   if(!diag)
#      d<-diag.remove(d)
   #Now, construct a model

   rm <- dim(d)[1]
#   allval <- array(data = 0,dim = c(3,n))
   allval <- array(data = 0,dim = c(rm,3,n))
   bnum <- array(data=0,dim=c(k))
   for(i in 1:rm){
      for (ac in 1:n) {
         bnum[b[ac]] = bnum[b[ac]]+1
         for(j in 1:n) {
            if(b[j]==b[ac]){
                allval[i,1,ac]=allval[i,1,ac]+d[i,ac,j]
#                allval[1,ac]=allval[1,ac]+d[i,ac,j]
            } else {
                allval[i,2,ac]=allval[i,2,ac]+d[i,ac,j]
#                allval[2,ac]=allval[2,ac]+d[i,ac,j]
            }
         }
      }
   }

   for (ac in 1:n) {
       allval[1,3,ac]=allval[1,1,ac]+n-bnum[b[ac]]-allval[1,2,ac]
   }

   block.membership<-b[ec$cluster$order]
   order.vector<-ec$cluster$order
   if(length(dim(allval))>2){
      blocked.data<-allval[1,,ec$cluster$order]
   }else{
      blocked.data<-allval[1,,ec$cluster$order]
   }
   if(dim(allval)[1]==1){
      allegiance.model<-allval[1,,]
      colnames(allegiance.model) = ec$glabels[ec$cluster$order]
      rownames(allegiance.model) = c("help", "hurt","allegiance")
   }else{
      allegiance.model<-allval
      dimnames(block.model)<-list(glabels,ec$glabels[ec$cluster$order],c("help", "hurt","allegiance"))
   }
   plabels<-plabels[ec$cluster$order]
 #  glabels<-glabels
 #  rlabels<-rlabels
   cluster.method<-ec$cluster.method
   equiv.fun<-ec$equiv.fun
   equiv.metric<-ec$metric

   o <- list( block.membership=block.membership, order.vector = order.vector, blocked.data=blocked.data,allegiance.model=allegiance.model,plabels=plabels,cluster.method=cluster.method,equiv.fun=equiv.fun,equiv.metric=equiv.metric,block.num=bnum)
   class(o)<-"allegiance"
   o
}





rshrinker<-function(X) {

   m <- nrow(X)
   labeldata <- array(0,c(m))
   group <- array(1,1)
   filler <- array(1,1)
   rownames(labeldata) <- rownames(X)
   dat <- numeric()
   dat <- rbind(X[1,])
   labeldata[1] = 1
   
   for(i in 2:m) {
       temp <- X[i,]
       rd = nrow(dat)
       match = 0
       for( j in 1:rd) {
          if( (sum(abs(temp-dat[j,])))==0) {
             labeldata[i] = j
             match = 1
             group[j] = group[j]+1
             j = rd+1
          } 
       }
       if(match ==0) {
             dat <- rbind(dat,temp)
             group <- rbind(group,filler)
             labeldata[i] = nrow(dat)
       }
   }
   rownames(dat) <- 1:dim(dat)[1]
   o<-list(data=dat,label=labeldata,group=group)
   class(o)<-"rshrunk"
   o
}
           
cshrinker<-function(X) {

   m <- ncol(X)
   labeldata <- array(0,c(1,m))
   group <- array(1,1)
   filler <- array(1,1)
   colnames(labeldata) <- colnames(X)
   dat <- numeric()
   dat <- cbind(X[,1])
   labeldata[1] = 1
   
   for(i in 2:m) {
       temp <- X[,i]
       rd = ncol(dat)
       match = 0
       for( j in 1:rd) {
          if( (sum(abs(temp-dat[,j])))==0) {
             labeldata[i] = j
             match = 1
             group[j] = group[j]+1
             j = rd+1
          } 
       }
       if(match ==0) {
             dat <- cbind(dat,temp)
             group <- cbind(group,filler)
             labeldata[i] = ncol(dat)
       }
   }
   colnames(dat) <- 1:dim(dat)[2]
   o<-list(data=dat,label=labeldata,group=group)
   class(o)<-"cshrunk"
   o
}

# CBIND
# temp = array(2,c(10))
# X = array(3,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# temp = array(3,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# temp = array(4,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# temp = array(3,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# temp = array(3,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# temp = array(4,c(10))
# X = cbind(X,temp)
# temp = array(2,c(10))
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)
# X = cbind(X,temp)



#####RBIND
# temp = array(2,c(10))
# X = array(3,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# temp = array(3,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# temp = array(4,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# temp = array(3,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# temp = array(3,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# temp = array(4,c(10))
# X = rbind(X,temp)
# temp = array(2,c(10))
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)
# X = rbind(X,temp)

data.row.rel <- function(X)
{
   # X is m by n with m actors and n events
   m <- ncol(X)
   for(j in 1:m)
      X[,j] <- X[,j]/sum(X[,j])
   # return thi
   X
}

data.col.rel <- function(X)
{
   # X is m by n with m actors and n events
   m <- nrow(X)
   for(j in 1:m)
      X[j,] <- X[j,]/sum(X[j,])
   # return thi
   X
}

autoprob <- function(X)
{
  X <- data.row.rel(X)
  X <- cshrinker(X)
  m <- nrow(X$data)
  n <- ncol(X$data)
  for(j in 1:n) {
     if(X$group[j] != 1)  {
        for(k in 1:m) {
           # binomial chance of no successes)
           X$data[k,j] = 1 - (1-X$data[k,j])^(X$group[j])
        }
     }
  }
  data <- X$data
  group <- X$group
  label <- X$label
  o <- list( data=data, group=group,label=label)
  o
}
             
coc<-function(A,k) {

  require(MASS)
  l = round(log(k,base=2))
  Done <- rowSums(A)
  Dtwo <- colSums(A)
  sgDone <- 1/sqrt(Done)
  sgDtwo <- 1/sqrt(Dtwo)
  C <- array(0,c(nrow(A),ncol(A)))
  An <- array(0,c(nrow(A),ncol(A)))
  for( q in 1:nrow(A))
     for( p in 1:ncol(A))
            C[q,p] = sgDone[q] * A[q,p]
  for( q in 1:nrow(A))
     for( p in 1:ncol(A))
         An[q,p] = C[q,p] * sgDtwo[p]
  Q = svd(An);
  Qu = as.matrix(Q$u[,2:(l+1)])
  Qv = as.matrix(Q$v[,2:(l+1)])
  dout = as.matrix(Qu)
  dtvt = as.matrix(Qv)
  rownames(dout) <- rownames(data)
  rownames(dtvt) <- colnames(data)
  for( q in 1:nrow(Qu))
     for( p in 1:ncol(Qu))
            dout[q,p] = sgDone[q] * Qu[q,p]
  for( q in 1:nrow(Qv))
     for( p in 1:ncol(Qv))
            dtvt[q,p] = sgDtwo[q] * Qv[q,p]
  ztwo = rbind(dout,dtvt)
  C  = kmeans(ztwo,k)
  D = cbind(ztwo,as.matrix(C$cluster))
  D
}
  
  
datafy    <- function(q)
{
m = nrow(q)
n = ncol(q)
z <- matrix(data=0,nrow=m+n,ncol=m+n)
source("comem.R")
z[1:n,1:n] <- com.col.rel(q)
z[(n+1):(m+n),(n+1):(m+n)] <- com.row.rel(q)
if(m>n) {
   for(i in 1:m)
       q[i,] <- (q[i,]-min(q[i,]))/(max(q[i,])-min(q[i,]))
} else {
   for(i in 1:n)
       q[,i] <- (q[,i]-min(q[,i]))/(max(q[,i])-min(q[,i]))
}
z[1:n,(n+1):(m+n)] <- t(q)
z[(n+1):(m+n),1:n] <- q
aaa <- dimnames(q)[[1]]
bbb <- dimnames(q)[[2]]
qqq <- c(aaa,bbb)
colnames(z) <- qqq
rownames(z) <- qqq
z
}

datafy2    <- function(q)
{
m = nrow(q)
n = ncol(q)
z <- matrix(data=0,nrow=m+n,ncol=m+n)
source("comem.R")
#z[1:n,1:n] <- com.col.rel(q)
#z[(n+1):(m+n),(n+1):(m+n)] <- com.row.rel(q)
#if(m>n) {
#   for(i in 1:m)
#       q[i,] <- (q[i,]-min(q[i,]))/(max(q[i,])-min(q[i,]))
#} else {
#   for(i in 1:n)
#       q[,i] <- (q[,i]-min(q[,i]))/(max(q[,i])-min(q[,i]))
#}
z[1:n,(n+1):(m+n)] <- t(q)
z[(n+1):(m+n),1:n] <- q
aaa <- dimnames(q)[[1]]
bbb <- dimnames(q)[[2]]
qqq <- c(aaa,bbb)
colnames(z) <- qqq
rownames(z) <- qqq
z
}