rdpgtst.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)

## Loading required package: igraph

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

require(clue)

## Loading required package: clue

require(raster)

## Loading required package: raster

## Loading required package: sp

require(foreach)

## Loading required package: foreach

## foreach: simple, scalable parallel programming from Revolution Analytics
## Use Revolution R for scalability, fault tolerance and more.
## http://www.revolutionanalytics.com

require(ggplot2)

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.4.0

require(VN)

## Loading required package: 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")
    
    frakn <- 300 # number of vertices in graph 1
    fraknp <- 300 # number of vertices in graph 2
    feat <- 20 # number of feature vectors
    rhovec <- seq(0,1,by=.1)
    svec <- c(1:9)
    srhov <- expand.grid(s=svec,rho=rhovec)
    lensr <- nrow(srhov)
    
    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 <- 100
    
    
    MK <- foreach(srl = 1:lensr)%:%
          foreach(le = 1:len2,.combine='rbind')%dopar%{
    
            s <- srhov[srl,]$s
            corr <- srhov[srl,]$rho
    
            setseed <- 12345*s + 19*le
      
            ## A randomly generated  graph
            g1 <- sample_dot_product(lpvs)
            A <- as.matrix(get.adjacency(g1,type="both"))
            
            B <- adjcorrH(A,Pmat,corr)
            g2 <- graph_from_adjacency_matrix(B)
    
            voi <- sample(c(1:(frakn)),1)
            Nhtmp <- unlist(ego(g1,h,nodes=voi,mindist=1)) # mindist=0: close, 1: open
            S <- sample(Nhtmp,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="rdpgsrhovaryNOplot.RData")
}else{
    load("rdpgsrhovaryNOplot.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(srhov,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=df1rstr, aes(x=s, y = MKnrkbar,colour=as.factor(rho))) + 
        scale_x_continuous(breaks = df1rstr$s) +  
        geom_line(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(plot.title = element_text(hjust = 0.5)) + 
        geom_errorbar(aes(ymin=X1, ymax=X2,aes=.7),width=0.5,size=0.75) + 
        theme(text=element_text(size=25)) + 
        ylab(expression(paste(tau,"(x')"))) + 
        xlab(expression(s[x])) +
        guides(color=guide_legend(title=expression(rho))) +
        labs(title= expression(paste("RDPG: Vary ", s[x], " and ", rho)))

## Warning: Ignoring unknown aesthetics: aes

#pdf("rdpg_srhovary.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