reports/multistart.R

   1 library(morpheus)
   2
   3 testMultistart <- function(N, n, K, p, beta, b, link, nstart, ncores)
   4 {
   5   ms <- multiRun(
   6     list(n=n,p=p,beta=beta,b=b,optargs=list(K=K,link=link,nstart=nstart)),
   7     list(
   8       function(fargs) {
   9         # 1 start
  10         library(morpheus)
  11         K <- fargs$optargs$K
  12         op <- optimParams(K, fargs$optargs$link, fargs$optargs)
  13         x_init <- list(p=rep(1/K,K-1), beta=fargs$mu, b=rep(0,K))
  14                 res <- NULL
  15                 tryCatch({
  16           res <- do.call(rbind, op$run(x_init))
  17                 }, error = function(e) {
  18                     res <- NA
  19                 })
  20                 res
  21       },
  22       function(fargs) {
  23         # B starts
  24         library(morpheus)
  25         K <- fargs$optargs$K
  26         op <- optimParams(K, fargs$optargs$link, fargs$optargs)
  27         best_val <- Inf
  28         best_par <- list()
  29         for (i in 1:fargs$optargs$nstart)
  30         {
  31           #x_init <- list(p=rep(1/K,K-1), beta=i*fargs$mu, b=rep(0,K))
  32           M <- matrix(rnorm(d*K), nrow=d, ncol=K)
  33           M <- t(t(M) / sqrt(colSums(M^2)))
  34           x_init <- list(p=rep(1/K,K-1), beta=M, b=rep(0,K))
  35           tryCatch({
  36             par <- op$run(x_init)
  37           }, error = function(e) {
  38             par <- NA
  39           })
  40           if (!is.na(par[0]))
  41           {
  42             val <- op$f( op$linArgs(par) )
  43             if (val < best_val)
  44             {
  45               best_par <- par
  46               best_val <- val
  47             }
  48           }
  49         }
  50         # Bet that at least one run succeded:
  51         do.call(rbind,best_par)
  52       }
  53     ),
  54     prepareArgs = function(fargs, index) {
  55       library(morpheus)
  56       io = generateSampleIO(fargs$n, fargs$p, fargs$beta, fargs$b, fargs$optargs$link)
  57       fargs$optargs$M <- computeMoments(io$X, io$Y)
  58       mu <- computeMu(io$X, io$Y, fargs$optargs)
  59       fargs$mu <- mu
  60       fargs
  61     }, N=N, ncores=ncores, verbose=TRUE)
  62   for (i in 1:2)
  63     ms[[i]] <- alignMatrices(ms[[i]], ref=rbind(p,beta,b), ls_mode="exact")
  64   ms
  65 }
  66
  67 #model = binomial
  68 K <- 2
  69 p <- .5
  70 b <- c(-.2, .5)
  71 # Default values:
  72 link = "logit"
  73 N <- 10
  74 d <- 2
  75 n <- 1e4
  76 ncores <- 1
  77 nstart <- 3 #nstart-1 random starting points for each MC run
  78
  79 cmd_args <- commandArgs()
  80 for (arg in cmd_args)
  81 {
  82     if (substr(arg,1,1)!='-')
  83     {
  84         spl <- strsplit(arg,'=')[[1]]
  85         if (spl[1] == "nc") {
  86             ncores <- as.integer(spl[2])
  87         } else if (spl[1] == "N") {
  88             N <- as.integer(spl[2])
  89         } else if (spl[1] == "n") {
  90             n <- as.integer(spl[2])
  91         } else if (spl[1] == "d") {
  92             d <- as.integer(spl[2])
  93         } else if (spl[1] == "link") {
  94             link <- spl[2]
  95         } else if (spl[1] == "nstart") {
  96             nstart <- spl[2]
  97         }
  98     }
  99 }
 100 betas <- list(
 101     matrix( c(1,-2, 3,1), ncol=K ), #d=2
 102     matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=K ), #d=5
 103     matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0), ncol=K ) ) #d=10
 104 beta <- betas[[ ifelse( d==2, 1, ifelse(d==5,2,3) ) ]]
 105
 106 ms <- testMultistart(N, n, K, p, beta, b, link, nstart, ncores)
 107 ms_params <- list("N"=N, "nc"=ncores, "n"=n, "K"=K, "d"=d, "link"=link,
 108     "p"=c(p,1-sum(p)), "beta"=beta, "b"=b, "nstart"=nstart)
 109
 110 save("ms", "ms_params", file="multistart.RData")