#! /usr/bin/Rscript
###
### Generating and plotting the induced subgraphs of the high school
### Facebook and Survey networks generated by the vertices the two
### graphs have in common (shared vertices).
###
### Heather Gaddy Patsolic
### 2017 <[email protected]>
### S.D.G
#
require(raster)
require(igraph)
require(clue)
require(foreach)
require(ggplot2)
require(reshape2)
require(VN)
data('HSgraphs')
## User: If false, loads rdata file. If true, runs full script.
run <- FALSE
if(run){
require(doMC); nCores <- 4
#require(doMPI)
# Register Cores for parallelism
if("package:doMPI" %in% search()){
cl <- startMPIcluster()
registerDoMPI(cl)} else {
registerDoMC(nCores)
}
hhop <- 1 #2
g <- 0.1
R <- 50
iter <- 20
ell <- 1
# for where this data was generated, see hscore_plots.r
############################################
set.seed(1234)
perm2 <- sample(vcount(HSfrgraphcore),vcount(HSfrgraphcore),replace=FALSE)
adj2a <- get.adjacency(HSfrgraphcore)
adj2 <- adj2a[perm2,perm2]
map <- cbind(1:vcount(HSfrgraphcore),perm2)
colnames(map) <- c("core1id","core2id")
map2 <- cbind(perm2,1:vcount(HSfrgraphcore))
sp <- sort(perm2,index.return=TRUE)
map2sort <- map2[sp$ix,]
g.frscrambled <- graph_from_adjacency_matrix(adj2,mode="max")
P <- t(diag(1:82)[perm2,])
r1 <- apply(P,1,function(x){which(x==1)})
###stricken r2 <- apply(P,2,function(x){which(x)==1})
#induced_subgraph(HSfrgraphcore, FRpermlist[,1],
# impl = c("auto"))
#all(adj2a==(t(diag(1,82)[perm2,])%*%adj2%*%diag(1,82)[perm2,])) ##TRUE
#SET <- foreach(vvi = 1:2)%dopar%{
SET <- foreach(vvi = 1:nrow(adj2))%dopar%{
# note, we can use vvi as voi because of identity
# map between two graphs of core vertices.
set.seed <- 537*vvi
voi <- vvi
# no need to call on labels because identity map
fb.voi <- map2sort[voi,1]
fr.voi <- map2sort[voi,2]
#fb.voi<- which(HSfbgr$label==voi)
#friends.voi <- which(HSfriendsgr$label==voi)
Nhtmp <- unlist(ego(HSfbgraphcore,hhop,nodes=voi,mindist=1)) # mindist=0: close, 1: open
imp2 <- 0
nrankxp0 <- 0
nrankxpBad <- 0
nrvois <- rep(NA,length(Nhtmp))
imp2_h1l2 <- 0
nrankxp0_h1l2 <- 0
nrankxpBad_h1l2 <- 0
nrvois_h1l2 <- rep(NA,length(Nhtmp))
for(si in 1:length(Nhtmp)){
seeds <- Nhtmp[si];
S <- matrix(map2sort[seeds,],ncol=2)
frtmp <- sum(S[,2]>fr.voi)+fr.voi
# ### BEGINDEBUG
#
# s1 <- unlist(ego(HSfbgraphcore,hhop,nodes=S[,1],mindist=1)) # mindist=0: close, 1: open
# s2c <- unlist(ego(HSfrgraphcore,hhop,nodes=seeds,mindist=1)) # mindist=0: close, 1: open
# s2r1 <- unlist(ego(g.frscrambled,hhop,nodes=r1[seeds],mindist=1)) # mindist=0: close, 1: open
# #r1[s2c]==s2r1
# ## stricken--s2r2 <- unlist(ego(g.frscrambled,hhop,nodes=r2[seeds],mindist=1)) # mindist=0: close, 1: open
# ### ENDDEBUG
L <- localnbd(fb.voi,S,HSfbgraphcore,g.frscrambled,hhop,ell,R,g)
if(L$case=="possible"){
if(frtmp %in% L$labelsGxp){
#if(length(L$Cxp)==1){
# rankvoi <- 0
#}else{
Ps <- L$P
Ps <- Ps[1:length(L$labelsGx),1:length(L$labelsGxp)]
rankPs <- pass2ranksuplus(Ps)
rankvoi <- rankPs[(nrow(S)+1),
which(L$labelsGxp==(frtmp))] #note, this assumes all seeds used which will be the case in our rigged example --- labeling issue comes in if not all seeds used in vnsgm algorithm
rankvoi <- (rankvoi - 1)/(length(L$Cxp)-1)
}else{
rankvoi <- NA
imp2 <- imp2 + 1
}
if(!is.na(rankvoi)){
if(rankvoi == 0){
nrankxp0 <- nrankxp0 + 1
}else{
if(rankvoi >= 0.5){
nrankxpBad <- nrankxpBad + 1
}
}
}
}else{
rankvoi <- NA
imp2 <- imp2 + 1
}
nrvois[si] <- rankvoi
L_h1l2 <- localnbd(fb.voi,S,HSfbgraphcore,g.frscrambled,hhop,(ell+1),R,g)
if(L_h1l2$case=="possible"){
if(frtmp %in% L_h1l2$labelsGxp){
#if(length(L_h1l2$Cxp)==1){
# rankvoi_h1l2 <- 0
#}else{
Ps <- L_h1l2$P
Ps <- Ps[1:length(L_h1l2$labelsGx),1:length(L_h1l2$labelsGxp)]
rankPs_h1l2 <- pass2ranksuplus(Ps)
rankvoi_h1l2 <- rankPs_h1l2[(nrow(S)+1),
which(L_h1l2$labelsGxp==(frtmp))] #note, this assumes all seeds used which will be the case in our rigged example --- labeling issue comes in if not all seeds used in vnsgm algorithm
rankvoi_h1l2 <- (rankvoi_h1l2 - 1)/(length(L_h1l2$Cxp)-1)
}else{
rankvoi_h1l2 <- NA
imp2_h1l2 <- imp2_h1l2 + 1
}
if(!is.na(rankvoi_h1l2)){
if(rankvoi_h1l2 == 0){
nrankxp0_h1l2 <- nrankxp0_h1l2 + 1
}else{
if(rankvoi_h1l2 >= 0.5){
nrankxpBad_h1l2 <- nrankxpBad_h1l2 + 1
}
}
}
}else{
rankvoi_h1l2 <- NA
imp2_h1l2 <- imp2_h1l2 + 1
}
nrvois_h1l2[si] <- rankvoi_h1l2
}
Nhtmp2 <- unlist(ego(HSfbgraphcore,(hhop+1),nodes=voi,mindist=1)) # mindist=0: close, 1: open
imp2_h2l2 <- 0
nrankxp0_h2l2 <- 0
nrankxpBad_h2l2 <- 0
nrvois_h2l2 <- rep(NA,length(Nhtmp2))
imp2_h2l3 <- 0
nrankxp0_h2l3 <- 0
nrankxpBad_h2l3 <- 0
nrvois_h2l3 <- rep(NA,length(Nhtmp2))
for(si2 in 1:length(Nhtmp2)){
seeds2 <- Nhtmp2[si2];
S2 <- matrix(map2sort[seeds2,],ncol=2)
frtmp <- sum(S2[,2]>fr.voi)+fr.voi
L_h2l2 <- localnbd(fb.voi,S2,HSfbgraphcore,g.frscrambled,(hhop+1),(ell+1),R,g)
if(L_h2l2$case=="possible"){
if(frtmp %in% L_h2l2$labelsGxp){
#if(length(L_h2l2$Cxp)==1){
# rankvoi_h2l2 <- 0
#}else{
Ps <- L_h2l2$P
Ps <- Ps[1:length(L_h2l2$labelsGx),1:length(L_h2l2$labelsGxp)]
rankPs_h2l2 <- pass2ranksuplus(Ps)
rankvoi_h2l2 <- rankPs_h2l2[(nrow(S2)+1),
which(L_h2l2$labelsGxp==(frtmp))] #note, this assumes all seeds used which will be the case in our rigged example --- labeling issue comes in if not all seeds used in vnsgm algorithm
rankvoi_h2l2 <- (rankvoi_h2l2 - 1)/(length(L_h2l2$Cxp)-1)
}else{
rankvoi_h2l2 <- NA
imp2_h2l2 <- imp2_h2l2 + 1
}
if(!is.na(rankvoi_h2l2)){
if(rankvoi_h2l2 == 0){
nrankxp0_h2l2 <- nrankxp0_h2l2 + 1
}else{
if(rankvoi_h2l2 >= 0.5){
nrankxpBad_h2l2 <- nrankxpBad_h2l2 + 1
}
}
}
}else{
rankvoi_h2l2 <- NA
imp2_h2l2 <- imp2_h2l2 + 1
}
nrvois_h2l2[si2] <- rankvoi_h2l2
L_h2l3 <- localnbd(voi,S2,HSfbgraphcore,HSfrgraphcore,(hhop+1),(ell+2),R,g)
if(L_h2l3$case=="possible"){
if(frtmp %in% L_h2l3$labelsGxp){
#if(length(L_h2l3$Cxp)==1){
# rankvoi_h2l3 <- 0
#}else{
Ps <- L_h2l3$P
Ps <- Ps[1:length(L_h2l3$labelsGx),1:length(L_h2l3$labelsGxp)]
rankPs_h2l3 <- pass2ranksuplus(Ps)
rankvoi_h2l3 <- rankPs_h2l3[(nrow(S2)+1),
which(L_h2l3$labelsGxp==(frtmp))] #note, this assumes all seeds used which will be the case in our rigged example --- labeling issue comes in if not all seeds used in vnsgm algorithm
rankvoi_h2l3 <- (rankvoi_h2l3 - 1)/(length(L_h2l3$Cxp)-1)
}else{
rankvoi_h2l3 <- NA
imp2_h2l3 <- imp2_h2l3 + 1
}
if(!is.na(rankvoi_h2l3)){
if(rankvoi_h2l3 == 0){
nrankxp0_h2l3 <- nrankxp0_h2l3 + 1
}else{
if(rankvoi_h2l3 >= 0.5){
nrankxpBad_h2l3 <- nrankxpBad_h2l3 + 1
}
}
}
}else{
rankvoi_h2l3 <- NA
imp2_h2l3 <- imp2_h2l3 + 1
}
nrvois_h2l3[si2] <- rankvoi_h2l3
}
list(nrankxp0 = nrankxp0, nrankxpBad = nrankxpBad,imp2=imp2,
potseeds=length(Nhtmp),numcand=length(L$Cxp),
Nh1=Nhtmp,Nh2=Nhtmp2,nrvois=nrvois,
nrvois_h1l2=nrvois_h1l2, imp2_h1l2=imp2_h1l2,
nrvois_h2l2=nrvois_h2l2, imp2_h2l2=imp2_h2l2,
nrvois_h2l3=nrvois_h2l3, imp2_h2l3=imp2_h2l3)
#)#,"s"=length(L$S))#,mstime)
}
#save.image("psummary_h12ell123prime.RData")
save.image("psummary_h12ell123.RData")
}else{
load("psummary_h12ell123.RData")
}
#########Generate Summary Plot
simp2 <- Reduce(c, lapply(SET,function(x){x$imp2_h2l2}))
snrankxp0 <- Reduce(c,lapply(SET,function(x){length(which(x$nrvois_h2l2==0))}))
snrankxpBad <- Reduce(c,lapply(SET,function(x){length(which(x$nrvois_h2l2>=0.5))}))
ssizehnbdx <- Reduce(c,lapply(SET,function(x){length(x$nrvois_h2l2)}))
#simp2 <- Reduce(c, lapply(SET,function(x){x$imp2}))
#snrankxp0 <- Reduce(c,lapply(SET,function(x){length(which(x$nrvois==0))}))
#snrankxpBad <- Reduce(c,lapply(SET,function(x){length(which(x$nrvois>=0.5))}))
#ssizehnbdx <- Reduce(c,lapply(SET,function(x){length(x$nrvois)}))
snrankbw <- ssizehnbdx - snrankxpBad - snrankxp0 - simp2
dfsum <- data.frame(voi = c(1:nrow(adj2)),snrankxp0=snrankxp0,
snrankbw = snrankbw,snrankxpBad=snrankxpBad,
simp2=simp2)
dfsum1 <- data.frame(voi = c(1:nrow(adj2)),"zero"=snrankxp0,
"between"=snrankbw,"chance"=snrankxpBad, "na"=simp2)
ds1 <- melt(dfsum1,id="voi")
colnames(ds1) <- c("voi","group","count")
ds1$groupf <- factor(ds1$group,levels=c("na","chance","between","zero"))
#levels(ds1$group) <- c("NA","chance","between","0")
cpPalpurple<- rev(c('#a6dba0','#008837','#c2a5cf','#7b3294'))
pAll <-
ggplot(ds1, aes(x=voi, y=count, group=groupf, fill=groupf)) +
scale_x_continuous(breaks=seq(0,80,by=5)) +
geom_bar(stat = 'identity') +
theme(text=element_text(size=20),legend.position="top") + #c(0.85,0.9)) +
#, legend.title=expression(tau)) +
#scale_fill_discrete(name =expression(tau)
scale_fill_manual(name =expression(paste(tau,"(x')")), values = cpPalpurple) +
ylab("count") +
xlab("x") +
labs(title="Summary of VN-SGM for each VOI")
#pdf("psummary_h2l2.pdf",width=10,height=6)
print(pAll)
#dev.off()
if(run){
if(exists("cl")){
closeCluster(cl)
mpi.quit()
}
}
#
## Time:
### Working status:
#### Comments:
#####Soli Deo Gloria