Created
May 5, 2021 09:29
-
-
Save aleszib/3c392e7380b99b1611288f55b6ba6793 to your computer and use it in GitHub Desktop.
R Code for linked stochastic blockmodeling
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(blockmodeling) | |
stochBlock<-function(M, clu, eps=0, each=FALSE,weights=NULL, eachProb=FALSE, nMode=NULL, limits=NULL, design0=TRUE, weightClusterSize = 1){ | |
wc<-weightClusterSize | |
n<-dim(M)[1] | |
# mean05<-function(x)mean(c(x,0.5), na.rm=TRUE) | |
if(is.null(weights)){ | |
weights<-M | |
weights[,]<-1 | |
} else if(any(dim(weights)!=dim(M))) stop("Weights have wrong dim!") | |
w<-weights | |
if(is.null(nMode)) nMode<-ifelse(is.list(clu),length(clu),1) | |
if(nMode>1){ | |
tmN<-sapply(clu,length) | |
clu<-lapply(clu,function(x)as.integer(factor(x))) | |
tmNclu<-sapply(clu,max) | |
for(iMode in 2:nMode){ | |
clu[[iMode ]]<-clu[[iMode ]]+sum(tmNclu[1:(iMode -1)]) | |
} | |
clu<-unlist(clu) | |
cluMode<-rep(1:nMode, times=tmN) | |
if(isTRUE(design0)){ | |
ignoreInd<-funByBlocks(M, clu=cluMode,FUN = function(x)var(as.vector(x), na.rm=TRUE))>0 | |
ignoreInd<-ignoreInd[cluMode, cluMode] | |
w<-w*ignoreInd | |
}else if(!is.null(design0)){ | |
if(dim(design0)==c(nMode,nMode)){ | |
ignoreInd<-design0[cluMode, cluMode] | |
w<-w*ignoreInd | |
} | |
} | |
Mode<-rep(1:nMode,times=tmNclu) | |
HBM<-funByBlocks(M, clu = cluMode,ignore.diag = TRUE, FUN = mean) | |
HBM<-HBM[Mode,Mode] | |
HBM[HBM==0]<-0.5 #to avoid NaN's latter on | |
} else{ | |
clu<-as.integer(factor(clu)) | |
tmNclu<-max(clu) | |
tmN<-length(clu) | |
HBM<-NULL | |
Mode<-rep(1:nMode,times=tmNclu) | |
} | |
funByBlocksAdj<-function (clu){ | |
dM <- dim(M) | |
nn <- ifelse(length(dM) == 2, 1, dM[1]) | |
k<-max(clu) ##assumes that clu includes numbers from 1:k | |
IM.V <- array(NA, dim = c(nn, k, k)) | |
# dimnames(IM.V) <- list(1:nn,1:k,1:k) | |
for (iNet in 1:nn) { | |
if (length(dM) == 3){ | |
iM <- M[iNet, , ] | |
iHBM<-HBM[iNet,,] | |
} else { | |
iM <- M | |
iHBM<-HBM | |
} | |
for (i in 1:k) { | |
for (j in 1:k) { | |
B <- iM[clu == i, clu == j, drop = FALSE] | |
if(i==j) diag(B) <- NA | |
IM.V[iNet, i, j] <- mean(c(B,iHBM[i,j]), na.rm=TRUE) | |
} | |
} | |
} | |
if (nn == 1) | |
return(IM.V[1, , ]) | |
else return(IM.V) | |
} | |
#debug(funByBlocksAdj) | |
llMax<-rep(-Inf, n) | |
clu<-as.integer(factor(clu)) | |
k<-max(clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
IM<-funByBlocksAdj(clu=clu) | |
tclu<-table(clu) | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
oldIM<-IM | |
oldIM[,]<-1 | |
MnoDiag<-M | |
diag(MnoDiag)<-NA | |
diag(w)<-0 #due to diagonal | |
#print(ll) | |
# off<-0.5/(n*n) | |
# IM<-ifelse(IM<off,off,IM) | |
# IM<-ifelse(IM>(1-off),(1-off),IM) | |
iMode<-1 | |
tresh<-tmN[1] | |
allowedClus<- 1:tmNclu[1] | |
best<-NULL | |
# best2<-NULL | |
#normalization of weights | |
w[w>0]<-w[w>0]/mean(w[w>0]) | |
ll<-wc*sum(tclu*lptclu)+sum(w*MnoDiag*log(IM[clu,clu]),na.rm=TRUE)+sum(w*(1-MnoDiag)*log(1-IM[clu,clu]),na.rm=TRUE) | |
oldLl<- -Inf | |
while(ll>oldLl){ | |
oldIM<-IM | |
oldClu<-clu | |
oldLl<-ll | |
for(i in 1:n){ | |
x<-c(M[i,-i],M[-i,i]) | |
iw<-c(w[i,-i],w[-i,i]) | |
#iw<-iw/mean(iw) #added to preserve the importance of lptclu[j] | |
lli<-rep(-Inf, k) | |
if(i>tresh){ | |
iMode<-iMode+1 | |
tresh<-tresh+tmN[iMode] | |
allowedClus<- (cumsum(tmNclu)[iMode-1]+1):cumsum(tmNclu)[iMode] | |
} | |
for(j in allowedClus){ | |
p<-c(IM[j,clu[-i]],IM[clu[-i],j]) | |
lli[j]<-wc*lptclu[j]+sum(iw*(x*log(p)+(1-x)*log(1-p))) | |
} | |
# cat(lli, "\n") | |
# cat("Old: ",clu[i],", New: ",which.max(lli),"\n", sep="") | |
llMax[i]<-max(lli) | |
#clu[i]<-which.max(lli) | |
#clu[i]<-sample(which(lli==llMax[i]),size = 1) | |
tmp<-which(lli==llMax[i]);if(length(tmp)>1)tmp<-sample(tmp,1) | |
# start diagnostics | |
# if(clu[i]!=tmp){ | |
# print(sprintf("i = %1$d, clu[i] = %2$d, tmp = %3$d", i, clu[i], tmp)) | |
# print(lli) | |
# } | |
# end diagnostics | |
clu[i]<-tmp | |
if(each){ | |
tclu<-table(clu) | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
if(length(tclu)<k){ | |
clu<-oldClu | |
IM<-oldIM | |
eps<-Inf | |
print("Too few clusters") | |
break | |
} | |
IM<-funByBlocksAdj(clu=clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
# IM<-ifelse(IM<off,off,IM) | |
# IM<-ifelse(IM>(1-off),(1-off),IM) | |
} | |
} | |
tclu<-table(clu) | |
### old version | |
# if(length(tclu)<k){ | |
# clu<-oldClu | |
# IM<-oldIM | |
# oldIM[,]<-1 | |
# if(each|!eachProb) break | |
# each<-TRUE | |
# tclu<-table(clu) | |
# ptclu<-tclu/tmN[Mode] | |
# lptclu<-log(ptclu) | |
# next | |
# } | |
while(length(tclu)<k){ | |
missingClus<-(1:k)[!(as.character(1:k) %in% names(table(clu)))] | |
iClu<-ifelse(length(missingClus)>1,sample(missingClus,size = 1),missingClus) | |
iMode<-Mode[iClu] | |
allowedClus<- (c(0,cumsum(tmNclu))[iMode]+1):cumsum(tmNclu)[iMode] | |
ids<-which(clu%in%allowedClus) | |
selectId<-ids[which.min(llMax[ids])] | |
clu[selectId]<-iClu | |
llMax[selectId]<-max(llMax) | |
tclu<-table(clu) | |
# diagnostics start | |
cat("Empty clusters:", missingClus,"\n") | |
cat("iClu:", iClu,"\n") | |
cat("Clusters table:\n", tclu,"\n") | |
# diagnostics end | |
} | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
IM<-funByBlocksAdj(clu=clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
# IM<-ifelse(IM<off,off,IM) | |
# IM<-ifelse(IM>(1-off),(1-off),IM) | |
ll<-wc*sum(tclu*lptclu)+sum(w*MnoDiag*log(IM[clu,clu]),na.rm=TRUE)+sum(w*(1-MnoDiag)*log(1-IM[clu,clu]),na.rm=TRUE) | |
# diagnostics start | |
# cat(sprintf("ll = %f, diff = %f, each = %s\n",ll, sum((oldIM-IM)^2),each)) | |
# diagnostics end | |
iMode<-1 | |
tresh<-tmN[1] | |
allowedClus<- 1:tmNclu[1] | |
if(ll<oldLl){ | |
if(is.null(best) || oldLl>best$ll){ | |
cat("Gooing wrong way, saveing best!\n") | |
cat("new best ll =",oldLl,"\n") | |
best<-list(ll=oldLl, IM=oldIM, clu=oldClu) | |
} | |
break ### This was added to stop optimization if ll stops improving | |
} | |
# if(sum((oldIM-IM)^2)<=eps) if(each|| !eachProb){ | |
# if(!is.null(best)){ | |
# # cat("Using old best!\n") | |
# # cat("restored ll =",best$ll,"\n") | |
# if(ll>best$ll){ | |
# best2<-list(clu=clu, IM=IM, ll=ll) | |
# } else best2<-FALSE | |
# clu<-best$clu | |
# IM<-best$IM | |
# ll<-best$ll | |
# best<-NULL | |
# each<-eachProb | |
# } else break | |
# } else each<-TRUE | |
} | |
#ll<-sum(tclu*lptclu)+sum(w*MnoDiag*log(IM[clu,clu]),na.rm=TRUE)+sum(w*(1-MnoDiag)*log(1-IM[clu,clu]),na.rm=TRUE) | |
# the code bellow (both very similar parts, once using best, once best2) makes sure that if somewhere in the process better than the last solution was found, it will be returned. This indicates problems with the search procedure! | |
if(!is.null(best)&&best$ll>ll){ | |
clu<-best$clu | |
IM<-best$IM | |
ll<-best$ll | |
# cat("Using best!\n") | |
} | |
# if((!is.null(best2)&&!identical(best2,FALSE))&&(best2$ll>ll)){ | |
# clu<-best2$clu | |
# IM<-best2$IM | |
# ll<-best2$ll | |
# # cat("Using best2!\n") | |
# } | |
# | |
cat(sprintf("ll = %f, weightClusterSize = %s\n",ll, wc)) | |
clu<-splitClu(clu=clu, n=tmN) | |
res<-list(M=M, clu=clu, IM=IM, n= tmN, err=-ll, ll=ll, best=list(list(M=M, clu=clu, IM=IM, n= tmN))) | |
# the best element is redundant and only included so that some functions (namely stochBlockORP) work without modifications. May be removed soon if these functions are updated. | |
class(res)<-"opt.par" | |
return(res) | |
} | |
stochBlockORP<-function(M, | |
k,#number of clusters/groups | |
rep,#number of repetitions/different starting partitions to check | |
save.initial.param=TRUE, #save the initial parametrs of this call | |
save.initial.param.opt=FALSE, #save the initial parametrs for calls to optParC or optParMultiC | |
deleteMs=TRUE, #delete networks/matrices from results of optParC or optParMultiC to save space | |
max.iden=10, #the maximum number of results that should be saved (in case there are more than max.iden results with minimal error, only the first max.iden will be saved) | |
switch.names=NULL,#should partitions that only differ in group names be considert equal (is c(1,1,2)==c(2,2,1)) | |
return.all=FALSE,#if 'FALSE', solution for only the best (one or more) partition/s is/are returned | |
return.ll=TRUE,#if 'FALSE', only the resoults of crit.fun are returned (a list of all (best) soulutions including errors), else the resoult is list | |
seed=NULL,#the seed for random generation of partitions | |
RandomSeed=NULL, # the state of .Random.seed (e.g. as saved previously). Should not be "typed" by the user | |
parGenFun = genRandomPar, #The function that will generate random partitions. It should accept argumetns: k (number of partitions by modes, n (number of units by modes), seed (seed value for random generation of partition), addParam (a list of additional parametres) | |
mingr=NULL, #minimal alowed group size (defaults to c(minUnitsRowCluster,minUnitsColCluster) if set, else to 1) - only used for parGenFun function | |
maxgr=NULL, #maximal alowed group size (default to c(maxUnitsRowCluster,maxUnitsColCluster) if set, else to Inf) - only used for parGenFun function | |
addParam=list( #list of additional parameters for gerenrating partitions. Here they are specified for dthe default function "genRandomPar" | |
genPajekPar = TRUE, #Should the partitions be generated as in Pajek (the other options is completly random) | |
probGenMech = NULL), #Here the probabilities for different mechanizems for specifying the partitions are set. If not set this is determined based on the previous parameter. | |
maxTriesToFindNewPar=rep*10, #The maximum number of partition try when trying to find a new partition to optimize that was not yet checked before | |
skip.par = NULL, #partitions to be skiped | |
printRep= ifelse(rep<=10,1,round(rep/10)), #should some information about each optimization be printed | |
n=NULL, #the number of units by "modes". It is used only for generating random partitions. It has to be set only if there are more than two modes or if there are two modes, but the matrix representing the network is onemode (both modes are in rows and columns) | |
nCores=1, #number of cores to be used 0 -means all available cores, can also be a cluster object, | |
useParLapply=FALSE, #should ply be used instead of foreach | |
cl = NULL, #the cluster to use (if formed beforehand) | |
stopcl = FALSE, # should the cluster be stoped | |
limits=NULL, #limits - a kind of pre-specified blockmodel | |
... #other paramters to stochBlock | |
){ | |
dots<-list(...) | |
if(is.null(switch.names)){ | |
switch.names<-is.null(dots$BLOCKS) | |
} | |
if(save.initial.param)initial.param<-c(tryCatch(lapply(as.list(sys.frame(sys.nframe())),eval),error=function(...)return("error")),dots=list(...))#saves the inital parameters | |
if(is.null(mingr)){ | |
if(is.null(dots$minUnitsRowCluster)){ | |
mingr<-1 | |
} else { | |
mingr<-c(dots$minUnitsRowCluster,dots$minUnitsColCluster) | |
} | |
} | |
if(is.null(maxgr)){ | |
if(is.null(dots$maxUnitsRowCluster)){ | |
maxgr<-Inf | |
} else { | |
maxgr<-c(dots$maxUnitsRowCluster,dots$maxUnitsColCluster) | |
} | |
} | |
nmode<-length(k) | |
res<-list(NULL) | |
ll<-NULL | |
dots<-list(...) | |
if(is.null(switch.names)){ | |
switch.names<-is.null(dots$BLOCKS) | |
} | |
if(save.initial.param)initial.param<-c(tryCatch(lapply(as.list(sys.frame(sys.nframe())),eval),error=function(...)return("error")),dots=list(...))#saves the inital parameters | |
if(is.null(n)) if(nmode==1){ | |
n<-dim(M)[1] | |
} else if(nmode==2){ | |
n<-dim(M)[1:2] | |
} else warning("Number of nodes by modes can not be determined. Parameter 'n' must be supplied!!!") | |
if(!is.null(RandomSeed)){ | |
.Random.seed <- RandomSeed | |
} else if(!is.null(seed))set.seed(seed) | |
on.exit({ | |
res1 <- res[which(ll==max(ll, na.rm = TRUE))] | |
best<-NULL | |
best.clu<-NULL | |
for(i in 1:length(res1)){ | |
for(j in 1:length(res1[[i]]$best)){ | |
if( | |
ifelse(is.null(best.clu), | |
TRUE, | |
if(nmode==1) ifelse(switch.names, | |
!any(sapply(best.clu,rand2,clu2=res1[[i]]$clu)==1), | |
!any(sapply(best.clu,function(x)all(x==res1[[i]]$clu))) | |
) else ifelse(switch.names, | |
!any(sapply(best.clu,function(x,clu2)rand2(unlist(x),clu2),clu2=unlist(res1[[i]]$clu))==1), | |
!any(sapply(best.clu,function(x)all(unlist(x)==unlist(res1[[i]]$clu)))) | |
) | |
) | |
){ | |
best<-c(best,res1[i]) | |
best.clu<-c(best.clu,list(res1[[i]]$clu)) | |
} | |
if(length(best)>=max.iden) { | |
warning("Only the first ",max.iden," solutions out of ",length(na.omit(ll))," solutions with maximal loglikelihood will be saved.\n") | |
break | |
} | |
} | |
} | |
names(best)<-paste("best",1:length(best),sep="") | |
if(any(na.omit(ll)==-Inf) || ss(na.omit(ll))!=0 || length(na.omit(ll))==1){ | |
cat("\n\nOptimization of all partitions completed\n") | |
cat(length(best),"solution(s) with maximal loglikelihood =", max(ll,na.rm=TRUE), "found.","\n") | |
}else { | |
cat("\n\nOptimization of all partitions completed\n") | |
cat("All",length(na.omit(ll)),"solutions have loglikelihood",ll[1],"\n") | |
} | |
call<-list(call=match.call()) | |
best<-list(best=best) | |
checked.par<-list(checked.par=skip.par) | |
if(return.all) res<-list(res=res) else res<-NULL | |
if(return.ll) ll<-list(ll=ll, err=-ll) else ll<-NULL | |
if(!exists("initial.param")){ | |
initial.param<-NULL | |
} else initial.param=list(initial.param) | |
res<-c(list(M=M),res,best,ll,checked.par,call,initial.param=initial.param, list(Random.seed=.Random.seed, cl=cl)) | |
class(res)<-"opt.more.par" | |
return(res) | |
}) | |
if(nCores==1||!require(doParallel)){ | |
if(nCores!=1) { | |
oldWarn<-options("warn") | |
options(warn=1) | |
warning("Only single core is used as package 'doParallel' or 'doRNG' (or both) is/are not available") | |
options(oldWarn) | |
} | |
for(i in 1:rep){ | |
if(printRep & (i%%printRep==0)) cat("\n\nStarting optimization of the partiton",i,"of",rep,"partitions.\n") | |
find.unique.par<-TRUE | |
ununiqueParTested=0 | |
while(find.unique.par){ | |
temppar<-parGenFun(n=n,k=k,mingr=mingr,maxgr=maxgr,addParam=addParam) | |
find.unique.par<- | |
ifelse(is.null(skip.par), | |
FALSE, | |
if(nmode==1) ifelse(switch.names, | |
any(sapply(skip.par,rand2,clu2=temppar)==1), | |
any(sapply(skip.par,function(x)all(x==temppar))) | |
) else ifelse(switch.names, | |
any(sapply(skip.par,function(x,clu2)rand2(unlist(x),clu2),clu2=unlist(temppar))==1), | |
any(sapply(skip.par,function(x)all(unlist(x)==unlist(temppar)))) | |
) | |
) | |
ununiqueParTested=ununiqueParTested+1 | |
endFun<-ununiqueParTested>=maxTriesToFindNewPar | |
if(endFun) { | |
break | |
} else if(ununiqueParTested%%10==0) cat(ununiqueParTested,"partitions tested for unique partition\n") | |
} | |
if(endFun) break | |
skip.par<-c(skip.par,list(temppar)) | |
if(printRep==1) cat("Starting partition:",blockmodeling:::unlistPar(temppar),"\n") | |
res[[i]]<-stochBlock(M=M, clu=temppar, ...) | |
if(deleteMs){ | |
res[[i]]$M<-NULL | |
} | |
ll[i]<-res[[i]]$ll | |
if(printRep==1) cat("Final loglikelihood:",ll[i],"\n") | |
if(printRep==1) cat("Final partition: ",blockmodeling:::unlistPar(res[[i]]$clu),"\n") | |
} | |
} else { | |
oneRep<-function(i,M,n,k,mingr,maxgr,addParam,rep,limits,...){ | |
if (printRep) cat("\n\nStarting optimization of the partiton", | |
i, "of", rep, "partitions.\n") | |
temppar<-parGenFun(n=n,k=k,mingr=mingr,maxgr=maxgr,addParam=addParam) | |
#skip.par<-c(skip.par,list(temppar)) | |
tres <- try(stochBlock(M=M, clu=temppar, ...)) | |
if(class(tres)=="try-error"){ | |
tres<-list("try-error"=tres, ll=-Inf, startPart=temppar) | |
} | |
if(deleteMs){ | |
tres$M<-NULL | |
} | |
return(list(tres)) | |
} | |
if(!require(doParallel)|!require(doRNG)) useParLapply<-TRUE | |
if(nCores==0){ | |
nCores<-detectCores()-1 | |
} | |
if(useParLapply) { | |
if(is.null(cl)) cl<-makeCluster(nCores) | |
clusterSetRNGStream(cl) | |
nC<-nCores | |
clusterExport(cl, varlist = c("stochBlock")) | |
clusterEvalQ(cl, expr=library(blockmodeling)) | |
res<-parLapplyLB(cl = cl,1:rep, fun = oneRep, M=M,n=n,k=k,mingr=mingr,maxgr=maxgr,addParam=addParam,rep=rep, limits=limits, ...) | |
if(stopcl) stopCluster(cl) | |
res<-lapply(res,function(x)x[[1]]) | |
} else { | |
library(doParallel) | |
library(doRNG) | |
if(is.null(cl)) cl<-makeCluster(nCores) | |
if(!getDoParRegistered()){ | |
registerDoParallel(cl) | |
} | |
nC<-getDoParWorkers() | |
res<-foreach(i=1:rep,.combine=c, .packages='blockmodeling', .export = "stochBlock") %dorng% oneRep(i=i,M=M,n=n,k=k,mingr=mingr,maxgr=maxgr,addParam=addParam,rep=rep, limits=limits,...) | |
if(stopcl) stopCluster(cl) | |
} | |
ll<-sapply(res,function(x)x$ll) | |
} | |
} | |
ll<-function(M, clu, weights=NULL, nMode=NULL, limits=NULL, design0=TRUE){ | |
n<-dim(M)[1] | |
# mean05<-function(x)mean(c(x,0.5), na.rm=TRUE) | |
if(is.null(weights)){ | |
weights<-M | |
weights[,]<-1 | |
} else if(any(dim(weights)!=dim(M))) stop("Weights have wrong dim!") | |
w<-weights | |
if(is.null(nMode)) nMode<-ifelse(is.list(clu),length(clu),1) | |
if(nMode>1){ | |
tmN<-sapply(clu,length) | |
clu<-lapply(clu,function(x)as.integer(factor(x))) | |
tmNclu<-sapply(clu,max) | |
for(iMode in 2:nMode){ | |
clu[[iMode ]]<-clu[[iMode ]]+sum(tmNclu[1:(iMode -1)]) | |
} | |
clu<-unlist(clu) | |
cluMode<-rep(1:nMode, times=tmN) | |
if(isTRUE(design0)){ | |
ignoreInd<-funByBlocks(M, clu=cluMode,FUN = function(x)var(as.vector(x), na.rm=TRUE))>0 | |
ignoreInd<-ignoreInd[cluMode, cluMode] | |
w<-w*ignoreInd | |
}else if(!is.null(design0)){ | |
if(dim(design0)==c(nMode,nMode)){ | |
ignoreInd<-design0[cluMode, cluMode] | |
w<-w*ignoreInd | |
} | |
} | |
Mode<-rep(1:nMode,times=tmNclu) | |
HBM<-funByBlocks(M, clu = cluMode,ignore.diag = TRUE, FUN = mean) | |
HBM<-HBM[Mode,Mode] | |
} else{ | |
clu<-as.integer(factor(clu)) | |
tmNclu<-max(clu) | |
tmN<-length(clu) | |
HBM<-NULL | |
Mode<-rep(1:nMode,times=tmNclu) | |
} | |
funByBlocksAdj<-function (clu){ | |
dM <- dim(M) | |
nn <- ifelse(length(dM) == 2, 1, dM[1]) | |
k<-max(clu) ##assumes that clu includes numbers from 1:k | |
IM.V <- array(NA, dim = c(nn, k, k)) | |
# dimnames(IM.V) <- list(1:nn,1:k,1:k) | |
for (iNet in 1:nn) { | |
if (length(dM) == 3){ | |
iM <- M[iNet, , ] | |
iHBM<-HBM[iNet,,] | |
} else { | |
iM <- M | |
iHBM<-HBM | |
} | |
for (i in 1:k) { | |
for (j in 1:k) { | |
B <- iM[clu == i, clu == j, drop = FALSE] | |
if(i==j) diag(B) <- NA | |
IM.V[iNet, i, j] <- mean(c(B,iHBM[i,j]), na.rm=TRUE) | |
} | |
} | |
} | |
if (nn == 1) | |
return(IM.V[1, , ]) | |
else return(IM.V) | |
} | |
llMax<-rep(-Inf, n) | |
clu<-as.integer(factor(clu)) | |
k<-max(clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
IM<-funByBlocksAdj(clu=clu) | |
tclu<-table(clu) | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
oldIM<-IM | |
oldIM[,]<-1 | |
MnoDiag<-M | |
diag(MnoDiag)<-NA | |
diag(w)<-0 #due to diagonal | |
#normalization of weights | |
w[w>0]<-w[w>0]/mean(w[w>0]) | |
llMat<-MnoDiag*log(IM[clu,clu])+(1-MnoDiag)*log(1-IM[clu,clu]) | |
ll<-sum(tclu*lptclu)+sum(w*llMat,na.rm=TRUE) | |
return(list(ll=ll, llMat=llMat, err=-ll)) | |
} | |
icl<-function(M, clu, type="bin"){ | |
# https://arxiv.org/pdf/1303.2962.pdf | |
if(class(M)%in% c("opt.more.par","opt.par")){ | |
clu<-blockmodeling::clu(M) | |
ll<-M$best$best1$ll | |
M<-M$M | |
} else { | |
if(missing(clu)) stop("Clu must be specified!") | |
ll<-ll(M,clu = clu) | |
} | |
if(is.list(clu)) stop("Only one-mode/level networks are implemeted") | |
K<-length(unique(clu)) | |
n<-nrow(M) | |
icl<-ll - 1/2*K^2*log(n*(n-1)) - (K-1)/2*log(n) | |
return(icl) | |
} | |
llStochBlockR<-function(M, clu, weights=NULL, nMode=NULL, limits=NULL, design0=TRUE, weightClusterSize = 1){ | |
wc<-weightClusterSize | |
n<-dim(M)[1] | |
if(is.null(weights)){ | |
weights<-M | |
weights[,]<-1 | |
} else if(any(dim(weights)!=dim(M))) stop("Weights have wrong dim!") | |
w<-weights | |
if(is.null(nMode)) nMode<-ifelse(is.list(clu),length(clu),1) | |
if(nMode>1){ | |
tmN<-sapply(clu,length) | |
clu<-lapply(clu,function(x)as.integer(factor(x))) | |
tmNclu<-sapply(clu,max) | |
for(iMode in 2:nMode){ | |
clu[[iMode ]]<-clu[[iMode ]]+sum(tmNclu[1:(iMode -1)]) | |
} | |
clu<-unlist(clu) | |
cluMode<-rep(1:nMode, times=tmN) | |
if(isTRUE(design0)){ | |
ignoreInd<-funByBlocks(M, clu=cluMode,FUN = function(x)var(as.vector(x), na.rm=TRUE))>0 | |
ignoreInd<-ignoreInd[cluMode, cluMode] | |
w<-w*ignoreInd | |
}else if(!is.null(design0)){ | |
if(dim(design0)==c(nMode,nMode)){ | |
ignoreInd<-design0[cluMode, cluMode] | |
w<-w*ignoreInd | |
} | |
} | |
Mode<-rep(1:nMode,times=tmNclu) | |
HBM<-funByBlocks(M, clu = cluMode,ignore.diag = TRUE, FUN = mean) | |
HBM<-HBM[Mode,Mode] | |
HBM[HBM==0]<-0.5 #to avoid NaN's latter on | |
} else{ | |
clu<-as.integer(factor(clu)) | |
tmNclu<-max(clu) | |
tmN<-length(clu) | |
HBM<-NULL | |
Mode<-rep(1:nMode,times=tmNclu) | |
} | |
funByBlocksAdj<-function (clu){ | |
dM <- dim(M) | |
nn <- ifelse(length(dM) == 2, 1, dM[1]) | |
k<-max(clu) ##assumes that clu includes numbers from 1:k | |
IM.V <- array(NA, dim = c(nn, k, k)) | |
# dimnames(IM.V) <- list(1:nn,1:k,1:k) | |
for (iNet in 1:nn) { | |
if (length(dM) == 3){ | |
iM <- M[iNet, , ] | |
iHBM<-HBM[iNet,,] | |
} else { | |
iM <- M | |
iHBM<-HBM | |
} | |
for (i in 1:k) { | |
for (j in 1:k) { | |
B <- iM[clu == i, clu == j, drop = FALSE] | |
if(i==j) diag(B) <- NA | |
IM.V[iNet, i, j] <- mean(c(B,iHBM[i,j]), na.rm=TRUE) | |
} | |
} | |
} | |
if (nn == 1) | |
return(IM.V[1, , ]) | |
else return(IM.V) | |
} | |
clu<-as.integer(factor(clu)) | |
k<-max(clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
IM<-funByBlocksAdj(clu=clu) | |
tclu<-table(clu) | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
MnoDiag<-M | |
diag(MnoDiag)<-NA | |
diag(w)<-0 #due to diagonal | |
#print(ll) | |
#normalization of weights | |
w[w>0]<-w[w>0]/mean(w[w>0]) | |
ll<-wc*sum(tclu*lptclu)+sum(w*MnoDiag*log(IM[clu,clu]),na.rm=TRUE)+sum(w*(1-MnoDiag)*log(1-IM[clu,clu]),na.rm=TRUE) | |
clu<-splitClu(clu=clu, n=tmN) | |
res<-list(M=M, clu=clu, IM=IM, n= tmN, err=-ll, ll=ll) | |
return(res) | |
} | |
llStochBlockRTest<-function(M, clu, weights=NULL, nMode=NULL, limits=NULL, design0=TRUE, weightClusterSize = 1, IM=NULL){ | |
wc<-weightClusterSize | |
n<-dim(M)[1] | |
if(is.null(weights)){ | |
weights<-M | |
weights[,]<-1 | |
} else if(any(dim(weights)!=dim(M))) stop("Weights have wrong dim!") | |
w<-weights | |
if(is.null(nMode)) nMode<-ifelse(is.list(clu),length(clu),1) | |
if(nMode>1){ | |
tmN<-sapply(clu,length) | |
clu<-lapply(clu,function(x)as.integer(factor(x))) | |
tmNclu<-sapply(clu,max) | |
for(iMode in 2:nMode){ | |
clu[[iMode ]]<-clu[[iMode ]]+sum(tmNclu[1:(iMode -1)]) | |
} | |
clu<-unlist(clu) | |
cluMode<-rep(1:nMode, times=tmN) | |
if(isTRUE(design0)){ | |
ignoreInd<-funByBlocks(M, clu=cluMode,FUN = function(x)var(as.vector(x), na.rm=TRUE))>0 | |
ignoreInd<-ignoreInd[cluMode, cluMode] | |
w<-w*ignoreInd | |
}else if(!is.null(design0)){ | |
if(dim(design0)==c(nMode,nMode)){ | |
ignoreInd<-design0[cluMode, cluMode] | |
w<-w*ignoreInd | |
} | |
} | |
Mode<-rep(1:nMode,times=tmNclu) | |
HBM<-funByBlocks(M, clu = cluMode,ignore.diag = TRUE, FUN = mean) | |
HBM<-HBM[Mode,Mode] | |
HBM[HBM==0]<-0.5 #to avoid NaN's latter on | |
} else{ | |
clu<-as.integer(factor(clu)) | |
tmNclu<-max(clu) | |
tmN<-length(clu) | |
HBM<-NULL | |
Mode<-rep(1:nMode,times=tmNclu) | |
} | |
funByBlocksAdj<-function (clu){ | |
dM <- dim(M) | |
nn <- ifelse(length(dM) == 2, 1, dM[1]) | |
k<-max(clu) ##assumes that clu includes numbers from 1:k | |
IM.V <- array(NA, dim = c(nn, k, k)) | |
# dimnames(IM.V) <- list(1:nn,1:k,1:k) | |
for (iNet in 1:nn) { | |
if (length(dM) == 3){ | |
iM <- M[iNet, , ] | |
iHBM<-HBM[iNet,,] | |
} else { | |
iM <- M | |
iHBM<-HBM | |
} | |
for (i in 1:k) { | |
for (j in 1:k) { | |
B <- iM[clu == i, clu == j, drop = FALSE] | |
if(i==j) diag(B) <- NA | |
IM.V[iNet, i, j] <- mean(c(B,iHBM[i,j]), na.rm=TRUE) | |
} | |
} | |
} | |
if (nn == 1) | |
return(IM.V[1, , ]) | |
else return(IM.V) | |
} | |
clu<-as.integer(factor(clu)) | |
k<-max(clu) | |
if(!is.null(limits)){ | |
for(i in 1:k)for(j in 1:k){ | |
IM[i,j] <- limits[[i,j]](IM[i,j]) | |
} | |
} | |
if(is.null(IM)) IM<-funByBlocksAdj(clu=clu) | |
tclu<-table(clu) | |
ptclu<-tclu/tmN[Mode] | |
lptclu<-log(ptclu) | |
MnoDiag<-M | |
diag(MnoDiag)<-NA | |
diag(w)<-0 #due to diagonal | |
#print(ll) | |
#normalization of weights | |
#w[w>0]<-w[w>0]/mean(w[w>0]) | |
w<-w/mean(w[w>0]) | |
ll<-wc*sum(tclu*lptclu)+sum(w*MnoDiag*log(IM[clu,clu]),na.rm=TRUE)+sum(w*(1-MnoDiag)*log(1-IM[clu,clu]),na.rm=TRUE) | |
clu<-splitClu(clu=clu, n=tmN) | |
res<-list(M=M, clu=clu, IM=IM, n= tmN, err=-ll, ll=ll, weights=w) | |
return(res) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment