rdpgratio.r

#! /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
###
### Heather Gaddy Patsolic 
### 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{
    load("rdpgratiovaryNOplot.RData")
}

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:
### Comments:
####Soli Deo Gloria