#! /usr/bin/Rscript
###
### Simulations for how the number of seeds, vois, and vertices
### in the graphs affects the accuracy of TopK.
### --
### Fix ratio of number of vertices in graphs: 1 (or .75)
### Fix number of vois: 1
### vary the number of seeds
###
### 2015 <[email protected]>
### S.D.G
#

require(igraph)
require(clue)
require(raster)
require(foreach)
require(ggplot2)
require(VN)

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)
}
#source("localnbd.r") # Was setting for plot, but this function is
#same as vnsgm.

frakn <- 300 # number of vertices in graph 1
fraknp <- 300 # number of vertices in graph 2
feat <- 20 # number of feature vectors
corr <- .6
svec <- 4
#svec <- c(2:5)
ratiovec <- seq(0.25,1,by=0.05)
rhovec <- 0.6
#rhovec <- c(0.6,0.7,0.8)
srv <- expand.grid(s=svec,rho=rhovec,ratio=ratiovec)
lensr <- nrow(srv)

set.seed(1234567)
# Random vectors for RDPG
lpvs <- matrix(rnorm(feat*frakn), feat, frakn) #  features, 10 vertices
lpvs <- apply(lpvs, 2, function(x) { return (abs(x)/sqrt(sum(x^2))) })
Pmat <- t(lpvs)%*%lpvs
#########################

v <- 1 # voi
h <- 2 # max walk from voi
ell <- 2
g <- 0.1 # gamma <- max tolerance for alpha
R <- 30
iternum <- 20
len2 <- 300

MK <- foreach(srl = 1:lensr)%:%
foreach(le = 1:len2,.combine='rbind')%dopar%{

s <- srv[srl,]$s r <- srv[srl,]$ratio
frakm <- frakn*r
corr <- srv[srl,]$rho set.seed <- 12345*srl + 19*le ## A randomly generated graph g1 <- sample_dot_product(lpvs) A <- as.matrix(get.adjacency(g1,type="both")) B <- adjcorrH(A,Pmat,corr) shared <- sample(c(1:frakn),frakm,replace=FALSE) unshared1 <- setdiff(c(1:frakn),shared) A <- A[c(shared,unshared1),c(shared,unshared1)] g1 <- graph_from_adjacency_matrix(A) B <- B[shared,shared] g2 <- graph_from_adjacency_matrix(B) voi <- sample(c(1:(frakm)),1) Nhtmp <- unlist(ego(g1,h,nodes=voi,mindist=1)) # mindist=0: close, 1: open S <- sample(intersect(Nhtmp,c(1:summ)),s) seedpot <- intersect(Nhtmp,c(1:frakm)) lsp <- length(seedpot) # what if lsp < s ??? S <- sample(seedpot,s) #L <- localnbd(voi,S,g1,g2,h,ell,R,g) seeds <- cbind(S,S) L <- vnsgm(voi,seeds,g1,g2,h,ell,R,g,sim=TRUE) if(L$case=="possible"){
Ps <- L$P Ps <- Ps[(1:length(L$labelsGx)),(1:length(L$labelsGxp))] rankPs <- pass2ranksuplus(Ps) rankvoi <- rankPs[(length(L$Sx)+1),(length(L$Sx)+1)] rankvoi <- (rankvoi - 1)/(length(L$Cxp)-1)
}else{
rankvoi <- NA
}

c(case=L$case,normrankvoi=rankvoi,"sx"=length(L$Sx),numcand=length(L$Cxp),"s"=length(L$S))#,mstime)

}

#save.image(file="rdpgratiovaryNOplot.RData")
}else{
}

MKnrkbar <- as.numeric(lapply(MK,function(x){mean(as.numeric(x[,2]))}))
MKnrksd <- as.numeric(lapply(MK,function(x){sd(as.numeric(x[,2]))}))
MKnrksd <- MKnrksd/sqrt(len2)
MKnrkci <- cbind(c(MKnrkbar-2*MKnrksd),c(MKnrkbar+2*MKnrksd))
# for 95% CI, use 1.96 instead of 2

df1 <- data.frame(srv,MKnrkbar, MKnrkci)

#rpref <- seq(0,1,by=0.2)
#rpref <- c(0,.3,.5,.7,1)
#rstr <- which(as.numeric(as.character(df1$r)) %in% rpref) #df1rstr <- df1[rstr,] p04 <- ggplot(data=df1, aes(x=ratio, y = MKnrkbar))+ #,colour=as.factor(rho))) + scale_x_continuous(breaks = df1$ratio) +
geom_line(size=1) + #aes(colour=as.factor(rho)),size=1)  +
geom_point() +
theme(axis.line.x = element_line(colour="black"),
axis.line.y = element_line(colour="black")) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_errorbar(aes(ymin=X1, ymax=X2,aes=.7),width=0.04,size=0.75) +
theme(text=element_text(size=25)) +
ylab(expression(paste(tau,"(x')"))) +
xlab(expression(r)) +
theme(plot.title = element_text(hjust = 0.5)) +
#guides(color=guide_legend(title="rho")) +
labs(title= "RDPG: Vary ratio (r)")
## Warning: Ignoring unknown aesthetics: aes
#        theme_bw() +
#        theme(plot.background = element_blank(),
#              panel.grid.major = element_blank(),
#              panel.grid.minor = element_blank()) +
#        theme(panel.border=element_blank()) +

#pdf("rdpg_ratiovary.pdf",width=8,height=6)
p04

#dev.off()

###save.image(file="srhovaryWITHplot.RData")

#if(run){
#    if(exists("cl")){
#        closeCluster(cl)
#        mpi.quit()
#    }
#}

#   Time:
##  Working status:
####Soli Deo Gloria