shvary.r

#! /usr/bin/Rscript
###
### Simulation for how h and ell
### affect the accuracy of VN-via-SGM.
### Fix number of vois: 1
###
### Heather Gaddy Patsolic
### 2017 <[email protected]>
### S.D.G
#

# CHECK doMC/doMPI, len2

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

require(ggplot2)

## Loading required package: ggplot2

require(VN)

## Loading required package: VN

## 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)
        }

    B = matrix(c(0.4,0.05,0.05,0.05,0.4,0.05,0.05,0.05,0.4),3,byrow=TRUE)
    n = c(100,100,100)

    np<-c(0,n)
    P <- matrix(0,sum(n),sum(n))
    for(i in 1:length(n)){
        for(j in 1:length(n)){
            P[(sum(np[1:i])+1):sum(np[1:i+1]),
              (sum(np[1:j])+1):sum(np[1:j+1])] <-
               # B[i,j]*matrix(1,np[i+1],np[j+1])
                B[i,j]*matrix(1,np[i+1],np[j+1])
        }
    }
    rho <- 0.6

    svec <- seq(1,9,by=1)

    m <- c(100,100,100)
    v <- 1 # number of voi
    hvec <- seq(1,4,by=1) # max walk from voi
    paird <- expand.grid(h=hvec,s=svec)
    lenprd <- nrow(paird)
    len2 <- 50
    sumn <- sum(n)

    MK <- foreach(hl = 1:lenprd)%:%
          foreach(le = 1:len2,.combine='rbind')%dopar%{
              set.seed(567*le+256*hl)

              s <- paird[hl,]$s
              S <- sample(c(1:sumn),s,replace=FALSE)
              h <- paird[hl,]$h
              voi <- sample(setdiff(c(1:sumn),S),v,replace=FALSE)

              setseed <- 456*le
              RM <- randmatcorr(B,n,n,P,rho,setseed)
              A <- RM$A
              g1 <- graph_from_adjacency_matrix(A)

              Nh <- unlist(ego(g1,h,nodes=voi,mindist=1)) # mindist=0: close, 1: open
              Sx1 <- sort(intersect(Nh,S)); sx <- length(Sx1)

              c(paird[hl,],"sx" = sx,"nh" = length(Nh))#,mstime)
    }
}else{
    load(file="handellvaryDataFrame.RData")
}

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

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

df1 <- data.frame(paird,sxbar, sxci)


p04 <-
        ggplot(data=df1, aes(x=s, y = sxbar),colour=as.factor(h)) +
        scale_x_continuous(breaks = df1$s) +
        scale_y_continuous(breaks = df1$s) +
        geom_line(aes(colour=as.factor(h)),size=1)  +
        geom_point(aes(colour=as.factor(h))) +
        theme(axis.line.x = element_line(colour="black"),
              axis.line.y = element_line(colour="black")) +
        geom_errorbar(aes(ymin=X1, ymax=X2,colour=as.factor(h)),
                      width=0.5,size=0.75) +
        theme(text=element_text(size=25)) +
        ylab(expression(s[x])) +
        xlab("s") +
        guides(color=guide_legend(title="h")) +
        labs(title= expression(paste(s[x]," as a function of h and s")))
p04

#pdf("shvary_sbm.pdf",width=8,height=8)
#p04
#dev.off()
#
#
#save.image(file="handellvaryDataFrame.RData")



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



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