R_temp/bipartite/LPA_wb_plus.R

In celiason/phenotools: Tools for the analysis and assembly of phenomic data in R

# LPA_wb_plus.R
# Label propagation algorithm for weighted bipartite networks that finds modularity.
# Contains the LPAwb+ and the DIRTLPAwb+ algorithms
# Author :  Stephen Beckett ( https://github.com/sjbeckett/weighted-modularity-LPAwbPLUS )
# MIT License

###############################################

LPA_wb_plus <- function(MATRIX,initialmoduleguess=NA) {
	
	#Make sure the smallest matrix dimension represent the red labels by making them the rows (If matrix is transposed here, will be transposed back at the end)
	flipped = 0
	if(dim(MATRIX)[1] > dim(MATRIX)[2]) {
		MATRIX = t(MATRIX)
		flipped = 1
	}

	Matsum = sum(MATRIX)
	col_marginals = colSums(MATRIX)
	row_marginals = rowSums(MATRIX)
	BMatrix = BarbersMatrix(MATRIX)

	#initiliase labels
	bluelabels = rep(NA,dim(MATRIX)[2])

	if (is.na(initialmoduleguess))
		redlabels = 1:dim(MATRIX)[1]
	else
		redlabels = sample(1:(initialmoduleguess+1),dim(MATRIX)[1],replace=TRUE)
	
	#Run Phase 1: Locally update lables to maximise Qb
	outlist = StageOne_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels)
	redlabels = outlist[[1]]
	bluelabels = outlist[[2]]
	Qb_now = outlist[[3]]

	#Run Phase 2: Connect divisions from top-down if it improves Qb, then run
	#phase 1 again. Repeat till Qb cannot be improved.
	outlist = StageTwo_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,Qb_now)
	redlabels = outlist[[1]]
	bluelabels = outlist[[2]]
	Qb_now = outlist[[3]]

	if(flipped==1) { #If matrix was flipped, swap row and column labels
    		holder = redlabels
		redlabels = bluelabels
		bluelabels = holder
	}

	return(list(Row_labels=redlabels, Col_labels=bluelabels, modularity=Qb_now))
}

###############################################

DIRT_LPA_wb_plus <- function(MATRIX,mini=4,reps=10) {
	A=LPA_wb_plus(MATRIX)

	mods=length(unique(A[[1]]))

	if((mods-mini) > 0) {
		for(aa in mini:mods) {
			for(bb in 1:reps) {
				B=LPA_wb_plus(MATRIX,aa)
				if(B[[3]]>A[[3]])
					A=B
			}
		}
	}

	return(list(Row_labels=A[[1]], Col_labels=A[[2]], modularity=A[[3]]))
}




###############################################

BarbersMatrix <- function(MATRIX) {
	return(MATRIX - (cbind(rowSums(MATRIX))%*%rbind(colSums(MATRIX)))/sum(MATRIX))
}

###############################################

WEIGHTEDMODULARITY <- function(BMatrix,Matsum,redlabels,bluelabels) {
	#see equation 8
	holdsum = 0

	for (rr in 1:length(redlabels)) {
	    for (cc in 1:length(bluelabels)) {
	        kroneckerdelta = redlabels[rr] == bluelabels[cc]
	        holdsum = holdsum + BMatrix[rr,cc] * kroneckerdelta
	    }
	}
	return(holdsum/Matsum)
}

###############################################
TRACE <- function(MATRIX) { return(sum(diag(MATRIX))) }

###############################################

WEIGHTEDMODULARITY2 <- function(BMatrix,Matsum,redlabels,bluelabels) {
	#see equation 9
	UNIred = unique(redlabels)
	Lred = length(UNIred)
	UNIblu = unique(bluelabels)
	Lblu = length(UNIblu)
	LABELMAT1 = matrix(0,Lred,length(redlabels))
	LABELMAT2 = matrix(0,length(bluelabels),Lblu)

	for(aa in 1:length(redlabels))
		LABELMAT1[which(UNIred == redlabels[aa]),aa] = 1

	for(aa in 1:length(bluelabels))
		LABELMAT2[aa,which(UNIblu == bluelabels[aa])] = 1

	return(TRACE(LABELMAT1 %*% BMatrix %*% LABELMAT2)/Matsum)

}

###############################################

DIVISION <- function(redlabels,bluelabels) {

	divisionsFound <- intersect(redlabels,bluelabels)

	return(divisionsFound)
}

###############################################

StageTwo_LPAwbdash <- function(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels, Qb_now) {
	
	divisionsFound = DIVISION(redlabels,bluelabels)
	NUMdiv = length(divisionsFound)
	IterateFlag = 1
	while(IterateFlag == 1) {
		CombinedDivisionsThisTime = 0
		if(NUMdiv > 1) {
			for(div1check in 1:(NUMdiv-1)) {
				Mod1 = divisionsFound[div1check]
				for(div2check in (div1check+1):NUMdiv) {
					CHECK_RED = redlabels
					CHECK_RED[which(redlabels==Mod1)] = divisionsFound[div2check]
					CHECK_BLUE = bluelabels
					CHECK_BLUE[which(bluelabels==Mod1)] = divisionsFound[div2check]


					QQ = WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECK_RED,CHECK_BLUE)			
					if(QQ > Qb_now) { #If a division will improve modularity - find the best way to do this
                    				FoundBetter = 0
                    				for(aa in 1:NUMdiv) {
                        				CHECK_RED2 = redlabels
                        				CHECK_RED2[which(redlabels==divisionsFound[aa])] = Mod1
                        				CHECK_BLUE2 = bluelabels
                        				CHECK_BLUE2[which(bluelabels==divisionsFound[aa])] = Mod1
                        				if(WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECK_RED2,CHECK_BLUE2) > QQ) {
                            					FoundBetter = 1
                        				}
                        				CHECK_RED2 = redlabels
                        				CHECK_RED2[which(redlabels==divisionsFound[aa])] = divisionsFound[div2check]
                        				CHECK_BLUE2 = bluelabels
                        				CHECK_BLUE2[which(bluelabels==divisionsFound[aa])] = divisionsFound[div2check]
                        				if(WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECK_RED2,CHECK_BLUE2) > QQ) {
                            					FoundBetter = 1
                        				}
                    				}
                    				if(FoundBetter == 0) { #If no better configuration found - JOIN.
                        				redlabels = CHECK_RED
                        				bluelabels = CHECK_BLUE
                        				CombinedDivisionsThisTime = CombinedDivisionsThisTime + 1
                    				}
                			}
				}
			}
			if(CombinedDivisionsThisTime == 0) {#If no divisions were joined move on
		            IterateFlag = 0
		        }
		}
		else {
		IterateFlag = 0		
		}

		outlist=StageOne_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels)  ##
		redlabels = outlist[[1]]
		bluelabels = outlist[[2]]
		Qb_now = outlist[[3]]
		divisionsFound = DIVISION(redlabels,bluelabels)
		NUMdiv = length(divisionsFound)
	}

return(list(redlabels, bluelabels, Qb_now))
}


###############################################

StageOne_LPAwbdash <- function(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels) {
	#Create storage containers for total marginals attached to each red(row)
	#label and blue(column) label

	BLUELABELLENGTH=length(bluelabels)
	REDLABELLENGTH=length(redlabels)
	TotalRedDegrees = rep(NA,max(redlabels))
	TotalBlueDegrees = rep(NA,max(BLUELABELLENGTH,REDLABELLENGTH))

	#Fill up these containers according to current labels
	#Red
	for(aa in 1:REDLABELLENGTH) {
	    if(is.na(TotalRedDegrees[redlabels[aa]])) {
        	TotalRedDegrees[redlabels[aa]] = row_marginals[aa]
	    }
	    else {
        	TotalRedDegrees[redlabels[aa]] = TotalRedDegrees[redlabels[aa]] + row_marginals[aa]
	    }
	}

	#Blue
	if(sum(is.na(bluelabels)) != BLUELABELLENGTH) { #occurs first time through as blue nodes unlabelled
	    for(bb in 1:BLUELABELLENGTH) {
        	if(is.na(TotalBlueDegrees[bluelabels[bb]])) {
	            TotalBlueDegrees[bluelabels[bb]] = col_marginals[bb]
		}
        	else {
        	    TotalBlueDegrees[bluelabels[bb]] = TotalBlueDegrees[bluelabels[bb]] + col_marginals[bb]
        	}
	    }
	}
	else {
	    TotalBlueDegrees = rep(0,max(BLUELABELLENGTH,REDLABELLENGTH))
	}
	
	
	#locally maximise modularity!
	outlist = LOCALMAXIMISATION(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,TotalRedDegrees,TotalBlueDegrees)
	redlabels = outlist[[1]]
	bluelabels = outlist[[2]]
	Qb_now = outlist[[3]]

	return(list(redlabels, bluelabels, Qb_now))

}



###############################################

LOCALMAXIMISATION <-  function(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,TotalRedDegrees,TotalBlueDegrees) {
	
	#Find score for current partition
	QbAfter = WEIGHTEDMODULARITY2(BMatrix,Matsum,redlabels,bluelabels)

	if(is.na(QbAfter)) { QbAfter = -999 }

	IterateFlag = 1
	while(IterateFlag == 1) {
		#Save old information
		QbBefore = QbAfter
		old_redlabels = redlabels
		old_bluelabels = bluelabels
		old_TRD = TotalRedDegrees
		old_TBD = TotalBlueDegrees

		#Update Blue Nodes using red node information (see equation 10)
		bluelabelchoices = unique(redlabels)
		
		for(bb in 1:length(bluelabels)) {
			if(is.na(bluelabels[bb]) == FALSE) {
				TotalBlueDegrees[bluelabels[bb]] = TotalBlueDegrees[bluelabels[bb]] - col_marginals[bb] 
			}
			changebluelabeltest = rep(NA,length(bluelabelchoices))

			for(ww in 1:length(bluelabelchoices)) {
				changebluelabeltest[ww] = sum( (redlabels == bluelabelchoices[ww]) * MATRIX[,bb])  -  col_marginals[bb]*TotalRedDegrees[bluelabelchoices[ww]]/Matsum 
			}

			#assign new label based on maximisation of above condition  

			labels = which(changebluelabeltest == max(changebluelabeltest,na.rm =TRUE))
			newlabelindex = labels[sample(1:length(labels),1)]
			bluelabels[bb] = bluelabelchoices[newlabelindex[1]]
			if(bluelabels[bb] > length(TotalBlueDegrees)) {
	                	TotalBlueDegrees[bluelabels[bb]] = 0
	        	}
            
        		#Update total marginals on new labelling
            		TotalBlueDegrees[bluelabels[bb]] = TotalBlueDegrees[bluelabels[bb]] + col_marginals[bb]
        	}

		#Now update red node labels based on blue node information (see equation 10)
		redlabelchoices = unique(bluelabels)
		

		for(aa in 1:length(redlabels)) {
			TotalRedDegrees[redlabels[aa]] = TotalRedDegrees[redlabels[aa]] - row_marginals[aa]
			changeredlabeltest = rep(NA,length(redlabelchoices))

			for(ww in 1:length(redlabelchoices)) {
				changeredlabeltest[ww] = sum( (bluelabels == redlabelchoices[ww]) * MATRIX[aa,])  -  row_marginals[aa]*TotalBlueDegrees[redlabelchoices[ww]]/Matsum  
			}
			
			#assign new label based on maximisation of above condition
			labels = which(changeredlabeltest == max(changeredlabeltest,na.rm = TRUE))
			newlabelindex = labels[sample(1:length(labels),1)]
 			redlabels[aa] = redlabelchoices[newlabelindex[1]]

			if(redlabels[aa] > length(TotalRedDegrees)) {
				TotalRedDegrees[redlabels[aa]] = 0
			}
			TotalRedDegrees[redlabels[aa]] = TotalRedDegrees[redlabels[aa]] + row_marginals[aa]
		}
		
		
		#Find the new modularity score based on node label updates.
        	QbAfter = WEIGHTEDMODULARITY(BMatrix,Matsum,redlabels,bluelabels)
        
        	#If this modularity is not as good as previous stop iterating and
        	#use that previous best information

        	if(QbAfter <= QbBefore) {
            		redlabels = old_redlabels
	            	bluelabels = old_bluelabels
        	    	TotalRedDegrees = old_TRD
        	    	TotalBlueDegrees = old_TBD
        	    	IterateFlag = 0
        	}
		
	}

	Qb_now = QbAfter
	
		
	return(list(redlabels, bluelabels, Qb_now))
}


###############################################