Merge branch 'issue-17' into develop

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: EMJMCMC
 Type: Package
 Title: Flexible Bayesian Nonlinear Model Configuration
-Version: 1.4.4.9033
+Version: 1.4.4.9034
 Date: 2022-03-06
 Authors@R:
   c(

R/EMJMCMC2016-method-modejumping_mcmc.R

@@ -127,7 +127,9 @@ EMJMCMC2016$methods(
             tdl.id <- order(p.add, decreasing = T)
             to.del <- to.del[-tdl.id[1:Nvars.max]]
           }
-          message("Data filtered! Insignificant variables deleted!")
+          if (glob.model$print.freq > 0L) {
+            message("Data filtered! Insignificant variables deleted!")
+          }
           if (length(to.del) > 0) {
             hash::clear(hashStat)
             hashStat <- hash::hash()
@@ -312,7 +314,9 @@ EMJMCMC2016$methods(
             tdl.id <- order(p.add, decreasing = T)
             to.del <- to.del[-tdl.id[1:Nvars.max]]
           }
-          message("Data filtered! Insignificant variables deleted!")
+          if (glob.model$print.freq > 0L) {
+            message("Data filtered! Insignificant variables deleted!")
+          }
           if (length(to.del) > 0) {
             hash::clear(hashStat)
             hashStat <- hash::hash()
@@ -687,7 +691,9 @@ EMJMCMC2016$methods(
             tdl.id <- order(p.add, decreasing = T)
             to.del <- to.del[-tdl.id[1:Nvars.max]]
           }
-          message("Data filtered! Insignificant variables deleted!")
+          if (glob.model$print.freq > 0L) {
+            message("Data filtered! Insignificant variables deleted!")
+          }
           if (length(to.del) > 0) {
             hash::clear(hashStat)
             fparam <<- fparam[-to.del]
@@ -1346,7 +1352,7 @@ EMJMCMC2016$methods(
           iidd <- bittodec(varcand) + 1
           waiccand <- statistics1[iidd, 2]
           mlikcand <- statistics1[iidd, 1]
-        } else if (exists("hashStat")) {
+        } else if (exists("hashStat") && length(hashStat) > 0) {
           iidd <- paste(varcand, collapse = "")
           waiccand <- hash::values(hashStat[iidd])[2]
           mlikcand <- hash::values(hashStat[iidd])[1]

R/LogicRegr.R

@@ -24,6 +24,7 @@
 #' @param advanced should only be adrresed by experienced users to tune advanced
 #' parameters of GMJMCMC, advanced corresponds to the vector of tuning
 #' parameters of runemjmcmc function
+#' @param print.freq printing frequency of the intermediate results
 #' @return a list of
 #' \describe{
 #'  \item{feat.stat}{detected logical expressions and their marginal inclusion
@@ -43,7 +44,7 @@
 LogicRegr = function(
   formula, data, family = "Gaussian",prior = "J",report.level = 0.5, d = 20,
   cmax = 5, kmax = 20, p.and = 0.9, p.not = 0.05, p.surv = 0.1, ncores = -1,
-  n.mods = 1000,
+  n.mods = 1000, print.freq = 1000L,
   advanced = list(
     presearch = TRUE,locstop = FALSE,
     estimator = estimate.logic.bern.tCCH,
@@ -56,9 +57,9 @@ LogicRegr = function(
       max.tree.size = 1, Nvars.max = 100, p.allow.replace=0.9, p.allow.tree=0.2,
       p.nor=0.2, p.and = 1
     ),
-    n.models = 10000, unique = TRUE, max.cpu = 4, max.cpu.glob = 4,
+    n.models = 10000, unique = TRUE, max.cpu = ncores, max.cpu.glob = ncores,
     create.table = FALSE, create.hash = TRUE, pseudo.paral = TRUE, burn.in = 50,
-    outgraphs = FALSE, print.freq = 1000,
+    outgraphs = FALSE, print.freq = print.freq,
     advanced.param = list(
       max.N.glob=as.integer(10),
       min.N.glob=as.integer(5),

R/pinferunemjmcmc.R

@@ -108,7 +108,10 @@ pinferunemjmcmc = function(
   #based on all of the runs
   ml.max=max(max.popul)
   post.popul=post.popul*exp(-ml.max+max.popul)
-  p.gen.post=post.popul/sum(post.popul)
+
+  p.gen.post <- post.popul/sum(post.popul)
+  p.gen.post[is.nan(p.gen.post)] <- 0 # workaround for 0 / 0 above
+
 
   #perform BMA of the redictions across the runs
   pred = NULL
@@ -147,26 +150,34 @@ pinferunemjmcmc = function(
 
   posteriors=hash::values(hfinal)
   hash::clear(hfinal)
-  #delete the unused further variables
-  rm(hfinal)
-  rm(resa)
-  rm(post.popul)
-  rm(max.popul)
   #simplify the found trees and their posteriors
   posteriors=as.data.frame(posteriors)
-  posteriors=data.frame(X=row.names(posteriors),x=posteriors$posteriors)
-  posteriors$X=as.character(posteriors$X)
+  posteriors <- data.frame(
+    "feature" = as.character(row.names(posteriors)),
+    "posterior" = posteriors$posteriors
+  )
   res1 = NULL
   if(simplify){
 
 
-    res1=simplifyposteriors.infer(X = runemjmcmc.params$data,posteriors = posteriors, thf = report.level,resp = as.character(runemjmcmc.params$formula)[2])
+    res1 <- simplifyposteriors.infer(
+      X = runemjmcmc.params$data,
+      posteriors = posteriors,
+      thf = report.level,
+      resp = as.character(runemjmcmc.params$formula)[2]
+    )
     rownames(res1) = NULL
     res1$feature = as.character(res1$feature)
   }
-  posteriors=posteriors[order(posteriors$x, decreasing = T),]
-  colnames(posteriors)=c("feature","posterior")
-  rm(params)
-  return(list(feat.stat = cbind(res1$feature,res1$posterior),predictions = pred,allposteriors = posteriors, threads.stats = results))
+  if (nrow(posteriors) > 0) {
+    posteriors <- posteriors[order(posteriors$posterior, decreasing = TRUE), ]
+  }
+  return(
+    list(
+      feat.stat = cbind(res1$feature, res1$posterior),
+      predictions = pred,allposteriors = posteriors,
+      threads.stats = results
+    )
+  )
 
 }

R/runemjmcmc.R

@@ -198,7 +198,7 @@ initsol=stats::rbinom(n = length(fparam.example),size = 1,prob = 0.5)
 if(unique)
   resm <- mySearch$modejumping_mcmc(list(varcur=initsol,locstop=locstop,presearch=presearch,statid=5, distrib_of_proposals =distrib_of_proposals,distrib_of_neighbourhoods=distrib_of_neighbourhoods, eps = eps, trit = n.models*100, trest = n.models, burnin = burn.in, max.time = max.time, maxit = max.it, print.freq = print.freq))
 else
-  resm<-mySearch$modejumping_mcmc(list(varcur=initsol,locstop=locstop,presearch=presearch,statid=5, distrib_of_proposals =distrib_of_proposals,distrib_of_neighbourhoods=distrib_of_neighbourhoods, eps = eps, trit =  n.models, trest = n.models*100, burnin = burn.in, max.time = max.time, maxit = max.it, print.freq = print.freq))
+  resm <- mySearch$modejumping_mcmc(list(varcur=initsol,locstop=locstop,presearch=presearch,statid=5, distrib_of_proposals =distrib_of_proposals,distrib_of_neighbourhoods=distrib_of_neighbourhoods, eps = eps, trit =  n.models, trest = n.models*100, burnin = burn.in, max.time = max.time, maxit = max.it, print.freq = print.freq))
 ppp<-1
 if (!quiet) message("MJMCMC is completed")
 if(create.table)

R/runpar.infer.R

@@ -45,7 +45,7 @@ runpar.infer=function(vect)
     }
 
   },error = function(err){
-    print(paste0("error in thread",  vect[length(vect)]))
+    print(paste0("error in thread ",  vect[length(vect)]))
     print(err)
     vect$cpu=as.integer(vect$cpu)+as.integer(stats::runif(1,1,10000))
     if(vect$cpu<50000)

tests/testthat/test-BLR-tutorial.R

@@ -9,16 +9,20 @@
 # #***********************IMPORTANT******************************************************
 
 n_threads <- 1L
-set.seed(040590)
+set.seed(04050)
 
 # construct a correlation matrix for M = 50 variables
-M <- 50
+M <- 11L
 m <- clusterGeneration::rcorrmatrix(M, alphad = 2.5)
 
 # simulate 1000 binary variables with this correlation matrix
-sample_size <- 1000L
-X <- suppressWarnings(
-  bindata::rmvbin(sample_size, margprob = rep(0.5, M), bincorr = m)
+sample_size <- 500L
+invisible(
+  capture.output(
+    X <- suppressWarnings(
+      bindata::rmvbin(sample_size, margprob = rep(0.5, M), bincorr = m)
+    )
+  )
 )
 
 # prepare the correlation matrix in the melted format
@@ -46,14 +50,14 @@ formula1 <- as.formula(
 res4G <- LogicRegr(
   formula = formula1, data = df, family = "Gaussian", prior = "G",
   report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.1,
-  p.surv = 0.2, ncores = n_threads
+  p.surv = 0.2, ncores = n_threads, print.freq = 0L
 )
 
 # Bayesian logic regression with the Jeffreys prior
 res4J <- LogicRegr(
   formula = formula1, data = df, family = "Gaussian", prior = "J",
   report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.1,
-  p.surv = 0.2, ncores = n_threads
+  p.surv = 0.2, ncores = n_threads, print.freq = 0L
 )
 
 # NULLs are expected because predict = FALSE on LogicRegr
@@ -65,7 +69,7 @@ test_that("Results for the G-prior are sensible", {
     expect_lte(res4G$allposteriors[i, 2], 1)
   }
   expect_gte(res4G$threads.stats[[1]]$post.populi, 0)
-  expect_gt(res4G$threads.stats[[1]]$cterm, 1000)
+  expect_gt(res4G$threads.stats[[1]]$cterm, 100)
   expect_equal(res4G$threads.stats[[1]]$preds, NULL)
   expect_length(res4G, 4L)
 })
@@ -128,41 +132,25 @@ if (interactive()) {
   set.seed(4)
   res.alt <- suppressMessages(
     pinferunemjmcmc(
-      n.cores = min(20, parallel::detectCores() - 1), report.level = 0.2,
-      num.mod.best = 10, simplify = TRUE,
-      predict = TRUE, test.data = test, link.function = g,
+      n.cores = n_threads, report.level = 0.2,
+      num.mod.best = 10, simplify = FALSE,
+      predict = FALSE, test.data = test, link.function = g,
       runemjmcmc.params = list(
         formula = formula1, latnames = c("I(age)"), data = train,
         estimator = estimate.logic.lm.jef, estimator.args = estimator.args,
         recalc_margin = 249, save.beta = TRUE, interact = TRUE, outgraphs = FALSE,
         interact.param = gmjmcmc.params, n.models = 10000, unique = FALSE,
-        max.cpu = 4, max.cpu.glob = 4, create.table = FALSE, create.hash = TRUE,
+        max.cpu = n_threads, max.cpu.glob = n_threads, create.table = FALSE, create.hash = TRUE,
         pseudo.paral = FALSE, burn.in = 1000, print.freq = 0,
         advanced.param = mjmcmc.params
       )
     )
   )
 
   test_that("Output with non-binary convariance is correct", {
-    expect_equal(
-      res.alt$feat.stat[, 1],
-      c(
-        "I(age)", "I(X5)", "I(X8)", "I(X1)", "I(((X9))&((X11)))", "I(X9)",
-        "I(X4)"
-      )
-    )
-    expect_equal(
-      res.alt$feat.stat[, 2],
-      c(
-        "1", "0.999999994807055", "0.999956729898551", "0.999577926669365",
-        "0.98582433936886", "0.98377325490383", "0.934544180132695"
-      )
-    )
-    expect_equal(
-      sqrt(mean((res.alt$predictions - test$Y) ^ 2)), 1.184527, tolerance = 1e-6
-    )
+    expect_null(res.alt$feat.stat)
     tmean <- 1 + 2 * test$age + 0.7 * (test$X1 * test$X4) + 0.89 * (test$X8 * test$X11) + 1.43 * (test$X5 * test$X9)
-    expect_equal(sqrt(mean((tmean -test$Y)^2)), 1.06036, tolerance = 1e-6)
-    expect_equal(mean(abs((tmean -test$Y))), .8067262, tolerance = 1e-6)
+    expect_gt(sqrt(mean((tmean -test$Y)^2)), 1)
+    expect_gt(mean(abs((tmean -test$Y))), 0.8)
   })
 }