R/the_mode_jumping_package2.0.r

Defines functions runemjmcmc simplify.formula estimate.glm.alt estimate.glm parallelize estimate.inla.ar1 estimate.inla.iid estimate.bas.lm.bagging pinferunemjmcmc mcgmjpse mcgmjpar LogicRegr estimate.logic.lm.jef estimate.logic.lm.tCCH estimate.logic.bern.tCCH estimate.logic.bern runpar.infer simplifyposteriors.infer estimate.logic.lm estimate.logic.lm.tCCH estimate.logic.bern.tCCH estimate.logic.bern estimate.logic.glm estimate.logic.lm estimate.bas.lm erf sigmoid estimate.bas.glm.bagging estimate.bigm estimate.speedglm estimate.elnet estimate.bas.glm.pen estimate.bas.glm

Documented in erf estimate.bas.glm estimate.bas.lm estimate.bigm estimate.elnet estimate.glm estimate.logic.glm estimate.logic.lm estimate.speedglm LogicRegr parallelize pinferunemjmcmc runemjmcmc sigmoid simplify.formula

#rm(list = ls(all = TRUE))

#install.packages("INLA", repos="http://www.math.ntnu.no/inla/R/testing")
#install.packages("bigmemory")
#install.packages("snow")
#install.packages("Rmpi")
#install.packages("ade4")
#install.packages("sp")
#install.packages("BAS")
#install.packages("hash")
#install.packages("stringi")
#install.packages("irlba")
#install.packages("bigalgebra")
#install.packages("speedglm")
#install.packages("biglm")
library(glmnet)
library(biglm)
library(hash)
library(sp)
library(INLA)
library(parallel)
library(bigmemory)
#library(snow)
library(MASS)
library(ade4)
#library(copula)
#library(compiler)
library(BAS)
library(stringi)
#library(speedglm)
require(stats)
#compile INLA


estimate.bas.glm <- function(formula, data, family, prior, logn)
{

  #only poisson and binomial families are currently adopted
  X <- model.matrix(object = formula,data = data)
  out <- bayesglm.fit(x = X, y = data[,1], family=family,coefprior=prior)
  # use dic and aic as bic and aic correspondinly
  return(list(mlik = out$logmarglik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + logn*out$rank),summary.fixed =list(mean = coefficients(out))))

}


estimate.bas.glm.pen <- function(formula, data, family, prior, logn,n,m,r=1)
{

  #only poisson and binomial families are currently adopted
  X <- model.matrix(object = formula,data = data)
  out <- bayesglm.fit(x = X, y = data[,1], family=family,coefprior=prior)
  fmla.proc<-as.character(formula)[2:3]
  fobserved <- fmla.proc[1]
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam <-stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj<-(stri_count_fixed(str = fparam, pattern = "&"))
  sj<-sj+(stri_count_fixed(str = fparam, pattern = "|"))
  Jprior <- sum(factorial(sj)/((m^sj)*2^(3*sj-2)))
  tn<-sum(stri_count_fixed(str = fmla.proc[2], pattern = "I("))

  return(list(mlik = (out$logmarglik+2*log(Jprior) + 2*tn*log(r)),waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + logn*out$rank),summary.fixed =list(mean = coefficients(out))))

}


#estimate elastic nets

estimate.elnet <- function(formula,response, data, family,alpha)
{
  X <- model.matrix(object = formula,data = data)
  if(dim(X)[2]<=2)
    return(list(mlik =  -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = array(0,dim=dim(X)[2]))))
  out <- glmnet(x=X[,-1], y = data[[response]], family=family, a=alpha)
  dout<- as.numeric(-deviance(out)[[out$dim[2]]])
  return(list(mlik = dout,waic = -dout , dic =  -dout,summary.fixed =list(mean = coef(out)[,out$dim[2]])))
}

# use dic and aic as bic and aic correspondinly
estimate.speedglm <- function(formula, data, family, prior) # weird behaviour, bad control of singularity
{
  X <- model.matrix(object = formula,data = data)
  out <- speedglm.wfit(y = data[,1], X = X, intercept=FALSE, family=family,eigendec = T, method = "Cholesky")
  if(prior == "AIC")
    return(list(mlik = -out$aic ,waic = -(out$deviance + 2*out$rank) , dic =  -(out$RSS),summary.fixed =list(mean = out$coefficients)))
  if(prior=="RSS")
    return(list(mlik = -out$RSS ,waic = -(out$deviance + 2*out$rank) , dic =  -(out$RSS),summary.fixed =list(mean = out$coefficients)))
}

estimate.bigm <- function(formula, data, family, prior, maxit = 2,chunksize = 1000000) # nice behaviour
{

  out <- bigglm(data = data, family=family,formula = formula, sandwich = F,maxit = maxit, chunksize = chunksize)
  if(prior == "AIC")
    return(list(mlik = -AIC(out,k = 2) ,waic = AIC(out,k = 2) , dic =  AIC(out,k = 2),summary.fixed =list(mean = coef(out))))
  if(prior=="RSS")
    return(list(mlik = -AIC(out,k = 2) ,waic = AIC(out,k = 2) , dic =  AIC(out,k = 2),summary.fixed =list(mean = coef(out))))
}

estimate.bas.glm.bagging <- function(formula, data, family, prior, logn, bag.size)
{

  #only poisson and binomial families are currently adopted
  X <- model.matrix(object = formula,data = data[sample.int(size = bag.size,n = dim(data)[1],replace = F),])
  out <- bayesglm.fit(x = X, y = data[,1], family=family,coefprior=prior)
  # use dic and aic as bic and aic correspondinly
  return(list(mlik = out$logmarglik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + logn*out$rank),summary.fixed =list(mean = coefficients(out))))

}

sigmoid<- function(x)
{
  return(1/(1+(exp(-x))))
}
erf <- function(x)
{
  return(2 * pnorm(x * sqrt(2)) - 1)
}
estimate.bas.lm <- function(formula, data, prior, n, g = 0)
{

  out <- lm(formula = formula,data = data)
  # 1 for aic, 2 bic prior, else g.prior

  p <- out$rank
  if(prior == 1)
  {
    ss<-sum(out$residuals^2)
    logmarglik <- -0.5*(log(ss)+2*p)
  }
  else if(prior ==2)
  {
    ss<-sum(out$residuals^2)
    logmarglik <- -0.5*(log(ss)+log(n)*p)
  }
  else
  {
    Rsquare <- summary(out)$r.squared
    #logmarglik =  .5*(log(1.0 + g) * (n - p -1)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
    logmarglik =  .5*(log(1.0 + g) * (n - p)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
  }

  # use dic and aic as bic and aic correspondinly
  return(list(mlik = logmarglik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coef(out))))

}

estimate.logic.lm <- function(formula, data, n, m, r = 1)
{
  out <- lm(formula = formula,data = data)
  p <- out$rank
  fmla.proc<-as.character(formula)[2:3]
  fobserved <- fmla.proc[1]
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam <-stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj<-(stri_count_fixed(str = fparam, pattern = "&"))
  sj<-sj+(stri_count_fixed(str = fparam, pattern = "|"))
  Jprior <- prod(factorial(sj)/((m^sj)*2^(3*sj-2)))
  #tn<-sum(stri_count_fixed(str = fmla.proc[2], pattern = "I("))
  mlik = (-BIC(out)+2*log(Jprior) + 2*p*log(r)+n)/2
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = mlik,waic = AIC(out)-n , dic =  BIC(out)-n,summary.fixed =list(mean = coef(out))))
}


estimate.logic.glm <- function(formula, data, family, n, m, r = 1)
{
  X <- model.matrix(object = formula,data = data)
  out <- bayesglm.fit(x = X, y = data[,1], family=family,coefprior=aic.prior())
  p <- out$rank
  fmla.proc<-as.character(formula)[2:3]
  fobserved <- fmla.proc[1]
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam <-stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj<-(stri_count_fixed(str = fparam, pattern = "&"))
  sj<-sj+(stri_count_fixed(str = fparam, pattern = "|"))
  Jprior <- sum(factorial(sj)/((m^sj)*2^(3*sj-2)))
  #tn<-sum(stri_count_fixed(str = fmla.proc[2], pattern = "I("))
  mlik = (out$logmarglik + 4*p - log(n)*p +2*log(Jprior) + 2*p*log(r)+n)
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = out$logmarglik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + log(n)*out$rank),summary.fixed =list(mean = coefficients(out))))
}
#define the function estimating parameters of a given Bernoulli logic regression with Jeffrey's prior
estimate.logic.bern = function(formula, data, family = binomial(), n=1000, m=50, r = 1,k.max=21)
{
  if(is.null(formula))
    return(list(mlik =  -10000 + rnorm(1,0,1),waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))

  out = glm(formula = formula,data = data, family=family)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = sum(log(factorial(sj)/((m^sj)*2^(2*sj-2))))
  mlik = (-(out$deviance + log(n)*(out$rank)) + 2*(Jprior))/2+n
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = mlik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + log(n)*out$rank),summary.fixed =list(mean = coefficients(out))))
}


#define the function estimating parameters of a given Bernoulli logic regression with robust g prior
estimate.logic.bern.tCCH = function(formula = NULL,y.id = 51, data, n=1000, m=50, r = 1, p.a = 1, p.b = 2, p.r = 1.5, p.s = 0, p.v=-1, p.k = 1,k.max=21)
{
  if(is.null(formula))
    return(list(mlik =  -10000 + rnorm(1,0,1),waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  X = scale(model.matrix(object = formula,data = data),center = T,scale = F)
  X[,1] = 1
  fmla.proc=as.character(formula)[2:3]
  out = glm(formula = as.formula(paste0(fmla.proc[1],"~X+0")),data=data,family = binomial())
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  beta=coef(out)[-1]
  if(length(which(is.na(beta)))>0)
  {
    return(list(mlik = -10000 + rnorm(1,0,1),waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  p.v = (n+1)/(p+1)
  sout = summary(out)
  J.a.hat = 1/sout$cov.unscaled[1,1]
  if(length(beta)>0&&length(beta)==(dim(sout$cov.unscaled)[1]-1)&&length(which(is.na(beta)))==0)
  {
    Q = t(beta)%*%solve(sout$cov.unscaled[-1,-1])%*%beta
  }else{
    return(list(mlik = -10000 + rnorm(1,0,1),waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  Jprior = sum(log(factorial(sj)/((m^sj)*2^(2*sj-2))))
  mlik = (logLik(out)- 0.5*log(J.a.hat) - 0.5*p*log(p.v) -0.5*Q/p.v + log(beta((p.a+p)/2,p.b/2)) + log(phi1(p.b/2,p.r,(p.a+p.b+p)/2,(p.s+Q)/2/p.v,1-p.k))+Jprior + p*log(r)+n)
  if(is.na(mlik)||mlik==-Inf)
    mlik = -10000+ rnorm(1,0,1)
  return(list(mlik = mlik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coefficients(out))))
}


#define the function estimating parameters of a given Gaussian logic regression with robust g prior
estimate.logic.lm.tCCH = function(formula = NULL, data, n=1000, m=50, r = 1, p.a = 1, p.b = 2, p.r = 1.5, p.s = 0, p.v=-1, p.k = 1,k.max=21)
{
  if(is.na(formula)||is.null(formula))
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  if(fmla.proc[2]=="-1")
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  out = lm(formula = formula,data = data)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = prod(factorial(sj)/((m^sj)*2^(2*sj-2)))
  p.v = (n+1)/(p+1)
  R.2 = summary(out)$r.squared

  mlik = (-0.5*p*log(p.v) -0.5*(n-1)*log(1-(1-1/p.v)*R.2) + log(beta((p.a+p)/2,p.b/2)) - log(beta(p.a/2,p.b/2)) + log(phi1(p.b/2,(n-1)/2,(p.a+p.b+p)/2,p.s/2/p.v,R.2/(p.v-(p.v-1)*R.2))) - hypergeometric1F1(p.b/2,(p.a+p.b)/2,p.s/2/p.v,log = T)+log(Jprior) + p*log(r)+n)
  if(mlik==-Inf||is.na(mlik)||is.nan(mlik))
    mlik = -10000
  return(list(mlik = mlik,waic = AIC(out)-n , dic =  BIC(out)-n,summary.fixed =list(mean = coef(out))))
}


#define the function estimating parameters of a given Gaussian logic regression with Jeffrey's prior
estimate.logic.lm = function(formula= NULL, data, n, m, r = 1,k.max=21)
{
  if(is.na(formula)||is.null(formula))
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  out = lm(formula = formula,data = data)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = prod(factorial(sj)/((m^sj)*2^(2*sj-2)))
  mlik = (-BIC(out)+2*log(Jprior) + 2*p*log(r)+n)/2
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = mlik,waic = AIC(out)-n , dic =  BIC(out)-n,summary.fixed =list(mean = coef(out))))
}


#define the function simplifying logical expressions at the end of the search
simplifyposteriors.infer=function(X,posteriors,th=0.0000001,thf=0.5, resp)
{
  todel = which(posteriors[,2]<th)
  if(length(todel)>0)
    posteriors=posteriors[-todel,]
  rhash=hash()
  for(i in 1:length(posteriors[,1]))
  {
    expr=posteriors[i,1]
    #print(expr)
    res=model.matrix(data=X,object = as.formula(paste0(resp,"~",expr)))
    res[,1]=res[,1]-res[,2]
    ress=c(stri_flatten(res[,1],collapse = ""),stri_flatten(res[,2],collapse = ""),posteriors[i,2],expr)
    ress[1] = stri_sub(ress[1],from = 1,to = 9999)
    if(!(ress[1] %in% values(rhash)||(ress[2] %in% values(rhash))))
      rhash[[ress[1]]]=ress
    else
    {
      if(ress[1] %in% keys(rhash))
      {
        rhash[[ress[1]]][3]= (as.numeric(rhash[[ress[1]]][3]) + as.numeric(ress[3]))
        if(stri_length(rhash[[ress[1]]][4])>stri_length(expr))
          rhash[[ress[1]]][4]=expr
      }
      else
      {
        rhash[[ress[2]]][3]= (as.numeric(rhash[[ress[2]]][3]) + as.numeric(ress[3]))
        if(stri_length(rhash[[ress[2]]][4])>stri_length(expr))
          rhash[[ress[2]]][4]=expr
      }
    }

  }
  res=as.data.frame(t(values(rhash)[c(3,4),]))
  res$V1=as.numeric(as.character(res$V1))
  tokeep = which(res$V1>thf)
  if(length(tokeep)>0)
  {
    res=res[tokeep,]
  }else
    warning(paste0("No features with posteriors above ",thf,". Returning everything"))
  res=res[order(res$V1, decreasing = T),]
  clear(rhash)
  rm(rhash)
  tokeep = which(res[,1]>1)
  if(length(tokeep)>0)
    res[tokeep,1]=1
  colnames(res)=c("posterior","feature")
  rownames(res) = NULL
  return(res)
}


#define a function performing the map step for a given thread
runpar.infer=function(vect)
{
  ret = NULL
  tryCatch({
    set.seed(as.integer(vect$cpu))
    do.call(runemjmcmc, vect[1:(length(vect)-5)])
    vals=values(hashStat)
    fparam=mySearch$fparam
    cterm=max(vals[1,],na.rm = T)
    ppp=mySearch$post_proceed_results_hash(hashStat = hashStat)
    post.populi=sum(exp(values(hashStat)[1,][1:vect$NM]-cterm),na.rm = T)

    betas = NULL
    mliks = NULL

    if(vect$save.beta){
      #get the modes of beta coefficients for the explored models
      Nvars=mySearch$Nvars
      linx =mySearch$Nvars+4
      lHash=length(hashStat)
      mliks = values(hashStat)[which((1:(lHash * linx)) %% linx == 1)]
      betas = values(hashStat)[which((1:(lHash * linx)) %% linx == 4)]
      for(i in 1:(Nvars-1))
      {
        betas=cbind(betas,values(hashStat)[which((1:(lHash * linx)) %% linx == (4+i))])
      }
      betas=cbind(betas,values(hashStat)[which((1:(lHash * linx)) %% linx == (0))])
    }

    preds = NULL
    if(vect$predict)
    {
      preds=mySearch$forecast.matrix.na(link.g = vect$link, covariates = (vect$test),betas = betas,mliks.in = mliks)$forecast
    }


    ret = list(post.populi = post.populi, p.post =  ppp$p.post, cterm = cterm, preds = preds, fparam = fparam, betas = betas, mliks = mliks )
    if(length(cterm)==0){
      vect$cpu=as.integer(vect$cpu)+as.integer(runif(1,1,10000))
      if(vect$cpu<50000)
        ret = runpar.infer(vect)
      else
        ret = NULL
    }

  },error = function(err){
    print(paste0("error in thread",  vect[length(vect)]))
    print(err)
    vect$cpu=as.integer(vect$cpu)+as.integer(runif(1,1,10000))
    if(vect$cpu<50000)
      ret = runpar.infer(vect)
    else
      ret =err
  },finally = {

    #clear(hashStat)
    #rm(hashStat)
    #rm(vals)
    gc()
    return(ret)

  })
}

#define the function estimating parameters of a given Bernoulli logic regression with Jeffrey's prior
estimate.logic.bern = function(formula, data, family = binomial(), n=1000, m=50, r = 1,k.max=21)
{
  if(is.null(formula))
    return(list(mlik =  -10000 + rnorm(1,0,1),waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))

  out = glm(formula = formula,data = data, family=family)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = sum(log(factorial(sj)/((m^sj)*2^(2*sj-2))))
  mlik = (-(out$deviance + log(n)*(out$rank)) + 2*(Jprior))/2+n
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = mlik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + log(n)*out$rank),summary.fixed =list(mean = coefficients(out))))
}


#define the function estimating parameters of a given Bernoulli logic regression with robust g prior
estimate.logic.bern.tCCH = function(formula = NULL,y.id = 51, data, n=1000, m=50, r = 1, p.a = 1, p.b = 2, p.r = 1.5, p.s = 0, p.v=-1, p.k = 1,k.max=21)
{
  if(is.null(formula))
    return(list(mlik =  -10000 + rnorm(1,0,1),waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  X = scale(model.matrix(object = formula,data = data),center = T,scale = F)
  X[,1] = 1
  fmla.proc=as.character(formula)[2:3]
  out = glm(formula = as.formula(paste0(fmla.proc[1],"~X+0")),data=data,family = binomial())
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  beta=coef(out)[-1]
  if(length(which(is.na(beta)))>0)
  {
    return(list(mlik = -10000 + rnorm(1,0,1),waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  p.v = (n+1)/(p+1)
  sout = summary(out)
  J.a.hat = 1/sout$cov.unscaled[1,1]
  if(length(beta)>0&&length(beta)==(dim(sout$cov.unscaled)[1]-1)&&length(which(is.na(beta)))==0)
  {
    Q = t(beta)%*%solve(sout$cov.unscaled[-1,-1])%*%beta
  }else{
    return(list(mlik = -10000 + rnorm(1,0,1),waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }

  Jprior = sum(log(factorial(sj)/((m^sj)*2^(2*sj-2))))
  mlik = (logLik(out)- 0.5*log(J.a.hat) - 0.5*p*log(p.v) -0.5*Q/p.v + log(beta((p.a+p)/2,p.b/2)) + log(phi1(p.b/2,p.r,(p.a+p.b+p)/2,(p.s+Q)/2/p.v,1-p.k))+Jprior + p*log(r)+n)
  if(is.na(mlik)||mlik==-Inf)
    mlik = -10000+ rnorm(1,0,1)
  return(list(mlik = mlik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coefficients(out))))
}


#define the function estimating parameters of a given Gaussian logic regression with robust g prior
estimate.logic.lm.tCCH = function(formula = NULL, data, n=1000, m=50, r = 1, p.a = 1, p.b = 2, p.r = 1.5, p.s = 0, p.v=-1, p.k = 1,k.max=21)
{
  if(is.na(formula)||is.null(formula))
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  if(fmla.proc[2]=="-1")
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  out = lm(formula = formula,data = data)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = prod(factorial(sj)/((m^sj)*2^(2*sj-2)))
  p.v = (n+1)/(p+1)
  R.2 = summary(out)$r.squared

  mlik = (-0.5*p*log(p.v) -0.5*(n-1)*log(1-(1-1/p.v)*R.2) + log(beta((p.a+p)/2,p.b/2)) - log(beta(p.a/2,p.b/2)) + log(phi1(p.b/2,(n-1)/2,(p.a+p.b+p)/2,p.s/2/p.v,R.2/(p.v-(p.v-1)*R.2))) - hypergeometric1F1(p.b/2,(p.a+p.b)/2,p.s/2/p.v,log = T)+log(Jprior) + p*log(r)+n)
  if(mlik==-Inf||is.na(mlik)||is.nan(mlik))
    mlik = -10000
  return(list(mlik = mlik,waic = AIC(out)-n , dic =  BIC(out)-n,summary.fixed =list(mean = coef(out))))
}


#define the function estimating parameters of a given Gaussian logic regression with Jeffrey's prior
estimate.logic.lm.jef = function(formula= NULL, data, n, m, r = 1,k.max=21)
{
  if(is.na(formula)||is.null(formula))
    return(list(mlik =  -10000,waic =10000 , dic =  10000,summary.fixed =list(mean = 1)))
  out = lm(formula = formula,data = data)
  p = out$rank
  if(p>k.max)
  {
    return(list(mlik = -10000,waic = 10000 , dic =  10000,summary.fixed =list(mean = 0)))
  }
  fmla.proc=as.character(formula)[2:3]
  fobserved = fmla.proc[1]
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]=stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam =stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]]
  sj=(stri_count_fixed(str = fparam, pattern = "&"))
  sj=sj+(stri_count_fixed(str = fparam, pattern = "|"))
  sj=sj+1
  Jprior = prod(factorial(sj)/((m^sj)*2^(2*sj-2)))
  mlik = (-BIC(out)+2*log(Jprior) + 2*p*log(r)+n)/2
  if(mlik==-Inf)
    mlik = -10000
  return(list(mlik = mlik,waic = AIC(out)-n , dic =  BIC(out)-n,summary.fixed =list(mean = coef(out))))
}


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 ,advanced = list(presearch = T,locstop = F ,estimator = estimate.logic.bern.tCCH,estimator.args =  list(data = data.example,n = 1000, m = 50,r=1),recalc_margin = 250, save.beta = F,interact = T,relations = c("","lgx2","cos","sigmoid","tanh","atan","erf"),relations.prob =c(0.4,0.0,0.0,0.0,0.0,0.0,0.0),interact.param=list(allow_offsprings=1,mutation_rate = 300,last.mutation = 5000, 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 = T,max.cpu = 4,max.cpu.glob = 4,create.table = F,create.hash = T,pseudo.paral = T,burn.in = 50,outgraphs=F,print.freq = 1000,advanced.param = list(
  max.N.glob=as.integer(10),
  min.N.glob=as.integer(5),
  max.N=as.integer(3),
  min.N=as.integer(1),
  printable = F)))
{
  advanced$formula = formula
  advanced$data = data
  advanced$interact.param$Nvars.max = d
  advanced$interact.param$ max.tree.size= cmax - 1
  advanced$interact.param$p.and = p.and
  advanced$interact.param$p.nor = p.not
  advanced$interact.param$p.allow.tree = p.surv
  if(!prior %in% c("J","G"))
  {
    warning("Wrong prior supplied. J (for Jeffrey's) and G (for robust g) priors are allowd only. Setting J as default.")
    prior = "J"
  }
  if(!family %in% c("Gaussian","Bernoulli"))
  {
    warning("Wrong familty supplied. Gaussian and Bernoulli families are allowd only. Setting Gaussian as default.")
    family = "Gaussian"
  }
  if(family == "Gaussian")
  {
    if(prior == "J")
    {
      advanced$estimator = estimate.logic.lm.jef
      advanced$estimator.args = list(data = data, n = dim(data)[1], m =stri_count_fixed(as.character(formula)[3],"+"),k.max = kmax)
    }else{
      advanced$estimator = estimate.logic.lm.tCCH
      advanced$estimator.args = list(data = data, n = dim(data)[1], m =stri_count_fixed(as.character(formula)[3],"+"),k.max = kmax)
    }
  }else{

    if(prior == "J")
    {
      advanced$estimator = estimate.logic.bern
      advanced$estimator.args =  list(data = data, n = dim(data)[1], m =stri_count_fixed(as.character(formula)[3],"+"),k.max = kmax)
    }else{
      advanced$estimator = estimate.logic.bern.tCCH
      advanced$estimator.args =  list(data = data, n = dim(data)[1], m =stri_count_fixed(as.character(formula)[3],"+"),k.max = kmax)
    }

  }

  if(ncores<1)
    ncores = detectCores()

  return(pinferunemjmcmc(n.cores = ncores,report.level =  report.level, simplify = T, num.mod.best = n.mods, predict = F, runemjmcmc.params = advanced))

}


mcgmjpar = function(X,FUN,mc.cores) mclapply(X= X,FUN = FUN,mc.preschedule = T,mc.cores = mc.cores)
mcgmjpse = function(X,FUN,mc.cores) lapply(X,FUN)

pinferunemjmcmc = function(n.cores = 4, mcgmj = mcgmjpar, report.level =  0.5, simplify = F, num.mod.best = 1000, predict = F,  test.data = 1, link.function = function(z)z, runemjmcmc.params)
{

  if(predict)
  {
    runemjmcmc.params$save.beta = T

    if(length(test.data)==0)
    {
      warning("Test data is not provided. No predictions will be made!")
    }
  }

  params = list(runemjmcmc.params)[rep(1,n.cores)]
  for(i in 1:n.cores)
  {
    params[[i]]$test = test.data
    params[[i]]$link = link.function
    params[[i]]$predict = predict
    params[[i]]$NM = num.mod.best
    params[[i]]$cpu=i
  }
  M = n.cores
  #results = runpar.infer(params[[1]])
  results=mcgmj(X = params,FUN = runpar.infer,mc.cores = n.cores)
  #print(results)
  #clean up
  gc()
  #prepare the data structures for final analysis of the runs
  compmax = runemjmcmc.params$interact.param$Nvars.max + 1
  resa=array(data = 0,dim = c(compmax,M*3))
  post.popul = array(0,M)
  max.popul = array(0,M)
  nulls=NULL
  not.null=1
  #check which threads had non-zero exit status
  for(k in 1:M)
  {
    if(length(results[[k]])<=1||length(results[[k]]$cterm)==0||length(results[[k]]$p.post)!=runemjmcmc.params$interact.param$Nvars.max)
    {
      nulls=c(nulls,k)
      #warning(paste0("Thread ",k,"did not converge or was killed by OS!"))
      next
    }
    else
    {
      not.null = k
    }

  }

  if(length(nulls) == M)
  {
    warning("All threads did not converge or gave an error! Returning stats from the threads only!")
    return(list(feat.stat = NULL,predictions = NULL,allposteriors = NULL, threads.stats = results))
  }


  #for all of the successful runs collect the results into the corresponding data structures
  for(k in 1:M)
  {
    if(k %in% nulls)
    {
      results[[k]]=results[[not.null]]
    }
    max.popul[k]=results[[k]]$cterm
    post.popul[k]=results[[k]]$post.populi
    resa[,k*3-2]=c(results[[k]]$fparam,"Post.Gen.Max")
    resa[,k*3-1]=c(results[[k]]$p.post,results[[k]]$cterm)
    resa[,k*3]=rep(post.popul[k],length(results[[k]]$p.post)+1)

  }
  #renormalize estimates of the marginal inclusion probabilities
  #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)

  #perform BMA of the redictions across the runs
  pred = NULL
  if(predict){
    pred = results[[1]]$preds*p.gen.post[1]
    for(i in 2:M)
    {

      pred=pred+results[[i]]$preds*p.gen.post[i]

    }
  }
  hfinal=hash()
  for(ii in 1:M)
  {
    resa[,ii*3]=p.gen.post[ii]*as.numeric(resa[,ii*3-1])
    resa[length(resa[,ii*3]),ii*3]=p.gen.post[ii]
    if(p.gen.post[ii]>0)
    {
      for(jj in 1:(length(resa[,ii*3])-1))
      {
        if(resa[jj,ii*3]>0)
        {
          if(as.integer(has.key(hash = hfinal,key =resa[jj,ii*3-2]))==0)
            hfinal[[resa[jj,ii*3-2]]]=as.numeric(resa[jj,ii*3])
          else
            hfinal[[resa[jj,ii*3-2]]]=hfinal[[resa[jj,ii*3-2]]]+as.numeric(resa[jj,ii*3])
        }

      }
    }
  }

  posteriors=values(hfinal)
  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)
  res1 = NULL
  if(simplify){


    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)
  gc()
  return(list(feat.stat = cbind(res1$feature,res1$posterior),predictions = pred,allposteriors = posteriors, threads.stats = results))

}

estimate.bas.lm.bagging <- function(formula, data, prior, n, g = 0, bag.size)
{

  out <- lm(formula = formula,data = data[sample.int(size = bag.size,n = dim(data)[1],replace = F),])
  # 1 for aic, 2 bic prior, else g.prior

  p <- out$rank
  if(prior == 1)
  {
    ss<-sum(out$residuals^2)
    logmarglik <- -0.5*(log(ss)+2*p)
  }
  else if(prior ==2)
  {
    ss<-sum(out$residuals^2)
    logmarglik <- -0.5*(log(ss)+log(n)*p)
  }
  else
  {
    Rsquare <- summary(out)$r.squared
    #logmarglik =  .5*(log(1.0 + g) * (n - p -1)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
    logmarglik =  .5*(log(1.0 + g) * (n - p)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
  }

  # use dic and aic as bic and aic correspondinly
  return(list(mlik = logmarglik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coef(out))))

}

estimate.inla.iid <- function(formula, args)
{

  out <- do.call(inla, c(args,formula = formula))
  # use dic and aic as bic and aic correspondinly
  coef<-out$summary.fixed$mode
  coef[1]<-coef[1]+out$summary.hyperpar$mode[1]
  return(list(mlik = out$logmarglik,waic = -(out$deviance + 2*out$rank) , dic =  -(out$deviance + logn*out$rank), summary.fixed =list(mean = coef)))

}

estimate.inla.ar1 <- function(formula, args)
{

  out<-NULL
  capture.output({withRestarts(tryCatch(capture.output({out <- do.call(inla, c(args,formula = formula)) })), abort = function(){onerr<-TRUE;out<-NULL})})
  if(is.null(out))
  {
    return(list(mlik = -10000,waic =10000, dic =10000, summary.fixed = 0))
  }
  # use dic and aic as bic and aic correspondinly
  coef<-out$summary.fixed$mode
  coef[1]<-coef[1]+out$summary.hyperpar$mode[1]/(1-out$summary.hyperpar$mode[2])
  return(list(mlik = out$mlik[1],waic = out$waic[1] , dic = out$dic[1], summary.fixed =list(mean =coef)))

}

parallelize<-function(X,FUN)
{
  max.cpu <- length(X)
  cl <-makeCluster(max.cpu ,type = paral.type,outfile = "")#outfile = ""
  clusterEvalQ(cl = cl,expr = c(library(INLA),library(bigmemory)))
  if(exists("statistics"))
  {
    clusterExport(cl=cl, "statistics")
    clusterEvalQ(cl,{statistics <- attach.resource(statistics);1})
  }
  res.par <- parLapply(cl = cl, X, FUN)
  stopCluster(cl)
  return(res.par)
}


estimate.glm <- function(formula, data, prior, family)
{

  out <- glm(formula = formula,data = data, family = family)
  # 1 for aic, 2 bic prior, else g.prior


  if(prior == 1)
  {

    logmarglik <- -AIC(out)
  }
  else
  {
    logmarglik <- -BIC(out)
  }

  # use dic and aic as bic and aic correspondinly
  return(list(mlik = logmarglik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coef(out))))

}


estimate.glm.alt <- function(formula, data, family, prior, n, g = 0)
{

  out <- glm(formula = formula, family = family, data = data)
  # 1 for aic, 2 bic prior, else g.prior

  p <- out$rank
  if(prior == 1)
  {

    logmarglik <- -AIC(out)
  }
  else if(prior ==2)
  {
    logmarglik <- -BIC(out)
  }
  else
  {
    Rsquare <- summary(out)$r.squared
    #logmarglik =  .5*(log(1.0 + g) * (n - p -1)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
    logmarglik =  .5*(log(1.0 + g) * (n - p)  - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
  }

  # use dic and aic as bic and aic correspondinly
  return(list(mlik = logmarglik,waic = AIC(out) , dic =  BIC(out),summary.fixed =list(mean = coef(out))))

}

simplify.formula<-function(fmla,names)
{
  fmla.proc<-as.character(fmla)[2:3]
  fobserved <- fmla.proc[1]
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = " ",replacement = "")
  fmla.proc[2]<-stri_replace_all(str = fmla.proc[2],fixed = "\n",replacement = "")
  fparam <- names[which(names %in% stri_split_fixed(str = fmla.proc[2],pattern = "+",omit_empty = F)[[1]] )]
  return(list(fparam = fparam,fobserved = fobserved))
}

# a function that creates an EMJMCMC2016 object with specified values of some parameters and deafault values of other parameters

runemjmcmc<-function(formula, data,
                     estimator,estimator.args = "list",n.models, unique = F,save.beta=F,latent="",max.cpu=4,max.cpu.glob=2,create.table=T, hash.length = 20, presearch=T, locstop =F ,pseudo.paral = F,interact = F,relations = c("","sin","cos","sigmoid","tanh","atan","erf"),relations.prob =c(0.4,0.1,0.1,0.1,0.1,0.1,0.1),interact.param=list(allow_offsprings=2,mutation_rate = 100,last.mutation=2000, max.tree.size = 10000, Nvars.max = 100, p.allow.replace = 0.7,p.allow.tree=0.1,p.nor=0.3,p.and = 0.7), recalc_margin = 2^10, create.hash=F,interact.order=1,burn.in=1, print.freq = 100,outgraphs=F,advanced.param=NULL, distrib_of_neighbourhoods=t(array(data = c(7.6651604,16.773326,14.541629,12.839445,2.964227,13.048343,7.165434,
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0.9936905,15.942490,11.040131,3.200394,15.349051,5.466632,14.676458,
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  1.5184551,9.285762,6.125034,3.627547,13.343413,2.923767,15.318774,
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  14.5295380,1.521960,11.804457,5.070282,6.934380,10.578945,12.455602,
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  6.0826035,2.453729,14.340435,14.863495,1.028312,12.685017,13.806295),dim = c(7,5))),  distrib_of_proposals = c(76.91870,71.25264,87.68184,60.55921,15812.39852))
{

  #first create the object
  assign("data.example",data, envir=globalenv())

  variables <- simplify.formula(formula,names(data.example))
  assign("fparam.example",variables$fparam, envir=globalenv())
  assign("fobserved.example",variables$fobserved, envir=globalenv())

  for(i in 1:length(fparam.example))
  {
    fparam.example[i]<<-paste("I(",variables$fparam[i],")",sep = "")
  }
  assign("mySearch",EMJMCMC2016(), envir=globalenv())
  mySearch$estimator <<- estimator
  mySearch$estimator.args <<- estimator.args
  mySearch$latent.formula <<- latent
  mySearch$save.beta <<- save.beta
  mySearch$recalc.margin <<- as.integer(recalc_margin)
  mySearch$max.cpu <<- as.integer(max.cpu)
  mySearch$locstop.nd <<- FALSE
  mySearch$max.cpu.glob <<- as.integer(max.cpu.glob)
  if(interact)
  {
    mySearch$allow_offsprings <<- as.integer(interact.param$allow_offsprings)
    mySearch$mutation_rate <<- as.integer(interact.param$mutation_rate)
    mySearch$Nvars.max <<- as.integer(interact.param$Nvars.max)
    mySearch$max.tree.size <<- as.integer(interact.param$max.tree.size)
    mySearch$p.allow.replace <<-  interact.param$p.allow.replace
    mySearch$p.allow.tree <<-  interact.param$p.allow.tree
    mySearch$sigmas<<-relations
    mySearch$sigmas.prob<<-relations.prob
    mySearch$p.nor <<- interact.param$p.nor
    mySearch$p.and <<- interact.param$p.and
    mySearch$last.mutation <<- as.integer(interact.param$last.mutation)
  }

  if(!is.null(advanced.param))
  {
    mySearch$max.N.glob<<-as.integer(advanced.param$max.N.glob)
    mySearch$min.N.glob<<-as.integer(advanced.param$min.N.glob)
    mySearch$max.N<<-as.integer(advanced.param$max.N)
    mySearch$min.N<<-as.integer(advanced.param$min.N)
    mySearch$printable.opt<<-advanced.param$printable
  }


  #distrib_of_proposals = с(0,0,0,0,10)
  if(exists("hashStat"))
  {
    clear(hashStat)
    remove(hashStat,envir=globalenv())
  }
  if(exists("statistics1"))
  {
    remove(statistics,envir=globalenv() )
    remove(statistics1,envir=globalenv())
  }
  if(exists("hash.keys1"))
  {
    remove(hash.keys,envir=globalenv())
    remove(hash.keys1,envir=globalenv())
  }

  if(create.table)
  {
    if(pseudo.paral)
      mySearch$parallelize <<- lapply
    #carry the search (training out)
    assign("statistics1",big.matrix(nrow = 2 ^min((length(fparam.example)),hash.length)+1, ncol =  16+length(fparam.example)*save.beta,init = NA, type = "double"), envir=globalenv())
    assign("statistics",describe(statistics1), envir=globalenv())

    mySearch$g.results[4,1]<<-0
    mySearch$g.results[4,2]<<-0
    mySearch$p.add <<- array(data = 0.5,dim = length(fparam.example))
    if((length(fparam.example))>20)
    {
      mySearch$hash.length<<-as.integer(hash.length)
      mySearch$double.hashing<<-T
      hash.keys1 <<- big.matrix(nrow = 2 ^(hash.length)+1, ncol = length(fparam.example),init = 0, type = "char")
      hash.keys <<- describe(hash.keys1)
    }

  }else if(create.hash)
  {

    assign("hashStat",hash(), envir=globalenv())
    mySearch$parallelize <<- lapply
    mySearch$hash.length<<-as.integer(20)
    mySearch$double.hashing<<-F
  }
  # now as the object is created run the algorithm
  initsol=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 = 0.000000001, trit = n.models*100, trest = n.models, burnin = burn.in, max.time = Inf, maxit = Inf, 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 = 0.000000001, trit =  n.models, trest = n.models*100, burnin = burn.in, max.time = Inf, maxit = Inf, print.freq = print.freq))
  ppp<-1
  print("MJMCMC is completed")
  if(create.table)
  {
    print("Post Proceed Results")
    ppp<-mySearch$post_proceed_results(statistics1 = statistics1)
    truth = ppp$p.post # make sure it is equal to Truth column from the article
    truth.m = ppp$m.post
    truth.prob = ppp$s.mass
    ordering = sort(ppp$p.post,index.return=T)
    print("pi truth")
    sprintf("%.10f",truth[ordering$ix])
    sprintf(fparam.example[ordering$ix])
  }
  else if(create.hash)
  {
    print("Post Proceed Results")
    ppp<-mySearch$post_proceed_results_hash(hashStat = hashStat)
    truth = ppp$p.post # make sure it is equal to Truth column from the article
    truth.m = ppp$m.post
    truth.prob = ppp$s.mass
    ordering = sort(ppp$p.post,index.return=T)
    print("pi truth")
    sprintf("%.10f",truth[ordering$ix])
    sprintf(fparam.example[ordering$ix])
  }

  if(outgraphs)
  {
    par(mar = c(10,4,4,2) + 4.1)
    barplot(resm$bayes.results$p.post,density = 46,border="black",main = "Marginal Inclusion (RM)",ylab="Probability",names.arg = mySearch$fparam,las=2)
    barplot(resm$p.post,density = 46,border="black",main = "Marginal Inclusion (MC)",ylab="Probability",names.arg = mySearch$fparam,las=2)
  }

  return(ppp)
}


# add plot(ppp), summary(ppp), print(ppp), ppp as a class itself, coef(ppp)

EMJMCMC2016 <- setRefClass(Class = "EMJMCMC2016",
                           fields = list(estimator.args = "list",
                                         max.cpu = "integer",
                                         objective = "integer",
                                         p.prior = "numeric",
                                         min.N = "integer",
                                         estimator = "function",
                                         parallelize = "function",
                                         parallelize.global = "function",
                                         parallelize.hyper  = "function",
                                         min.N.randomize = "integer",
                                         max.N.randomize = "integer",
                                         type.randomize = "integer",
                                         thin_rate = "integer",
                                         aa = "numeric",
                                         cc = "numeric",
                                         printable.opt = "logical",
                                         fobserved = "vector",
                                         switch.type = "integer",
                                         n.size ="integer",
                                         LocImprove = "array",
                                         max.N = "integer",
                                         save.beta = "logical",
                                         recalc.margin = "numeric",
                                         max.N.glob = "integer",
                                         min.N.glob = "integer",
                                         max.cpu.glob = "integer",
                                         max.cpu.hyper = "integer",
                                         switch.type.glob = "integer",
                                         isobsbinary = "array",
                                         fparam = "vector",
                                         fparam.pool = "vector",
                                         p.add = "array",
                                         latent.formula = "character",
                                         Nvars = "integer",
                                         seed = "integer",
                                         M.nd = "integer",
                                         locstop.nd = "logical",
                                         M.mcmc = "integer",
                                         SA.param = "list",
                                         p.allow.tree = "numeric",
                                         p.nor = "numeric",
                                         p.and = "numeric",
                                         sigmas.prob ="numeric",
                                         sigmas = "vector",
                                         p.allow.replace = "numeric",
                                         max.tree.size = "integer",
                                         double.hashing = "logical",
                                         hash.length = "integer",
                                         update.marg.mc = "logical",
                                         Nvars.max = "integer",
                                         Nvars.init = "integer",
                                         last.mutation = "integer",
                                         allow_offsprings = "integer",
                                         mutation_rate = "integer",
                                         g.results = "big.matrix"
                           ),
                           methods = list(
                             #class constructor
                             initialize = function(estimator.function = inla, estimator.args.list = list(family = "gaussian",data = data.example, control.fixed=list(prec=list(default= 0.00001),prec.intercept = 0.00001,mean=list(default= 0),mean.intercept = 0),
                                                                                                         control.family = list(hyper = list(prec = list(prior = "loggamma",param = c(0.00001,0.00001),initial = 0))),
                                                                                                         control.compute = list(dic = TRUE, waic = TRUE, mlik = TRUE)), search.args.list = NULL,latent.formula = "")
                             {
                               estimator <<- estimator.function
                               estimator.args <<- estimator.args.list
                               latent.formula <<- latent.formula
                               g.results <<- big.matrix(nrow = 4,ncol = 2)
                               g.results[1,1]<- -Inf
                               g.results[1,2]<- 1
                               g.results[2,1]<- Inf
                               g.results[2,2]<- 1
                               g.results[3,1]<- Inf
                               g.results[3,2]<- 1
                               g.results[4,1]<- 0
                               g.results[4,2]<- 0

                               if(is.null(search.args.list))
                               {
                                 max.cpu <<- as.integer(Nvars*0.05 + 1)
                                 objective <<- as.integer(1)
                                 if(Sys.info()['sysname']=="Windows")
                                 {
                                   parallelize <<- lapply
                                   parallelize.global <<- lapply
                                   parallelize.hyper  <<- lapply
                                 }
                                 else
                                 {
                                   parallelize <<- mclapply
                                   parallelize.global <<- mclapply
                                   parallelize.hyper  <<- mclapply
                                 }
                                 Nvars  <<- as.integer(length(fparam.example))
                                 min.N <<- as.integer(Nvars/6)
                                 min.N.glob <<- as.integer(Nvars/3)
                                 max.N.glob <<- as.integer(Nvars/2)
                                 max.N <<- as.integer(Nvars/5)
                                 switch.type.glob <<- as.integer(2)
                                 min.N.randomize <<- as.integer(1)
                                 max.N.randomize <<- as.integer(1)
                                 type.randomize <<- as.integer(3)
                                 max.cpu.glob <<- as.integer(Nvars*0.05 + 1)
                                 max.cpu.hyper <<- as.integer(2)
                                 save.beta <<- FALSE
                                 printable.opt <<- FALSE
                                 thin_rate<<- as.integer(-1)
                                 p.allow.tree <<- 0.6
                                 p.allow.replace <<- 0.3
                                 sigmas<<-c("","sin","cos","sigmoid","tanh","atan","erf")
                                 sigmas.prob<<-c(0.4,0.1,0.1,0.1,0.1,0.1,0.1)
                                 p.nor <<- 0.3
                                 p.and <<- 0.7
                                 max.tree.size<<- as.integer(15)
                                 Nvars.max <<- as.integer(Nvars)
                                 Nvars.init <<- as.integer(Nvars)
                                 allow_offsprings <<- as.integer(0)
                                 mutation_rate <<- as.integer(100)
                                 locstop.nd <<-FALSE
                                 double.hashing <<- (Nvars > 20)
                                 hash.length <<- as.integer(25)
                                 aa <<- 0.9
                                 cc <<- 0.0
                                 M.nd <<- as.integer(Nvars)
                                 M.mcmc <<- as.integer(5)
                                 SA.param <<- list(t.min = 0.0001,t.init = 10, dt = 3, M = as.integer(Nvars/5+1))
                                 fobserved <<- fobserved.example
                                 switch.type <<- as.integer(2)
                                 n.size <<- as.integer(10)
                                 LocImprove <<- as.array(c(50,50,50,50,150))
                                 isobsbinary <<- as.array(0:(length(fparam.example)-1))
                                 fparam <<- fparam.example
                                 fparam.pool <<- fparam.example
                                 p.add <<- array(data = 0.5,dim = Nvars)
                                 if(exists("statistics"))
                                 {
                                   recalc.margin <<- 2^Nvars
                                 }else if(exists("statistics1"))
                                 {
                                   recalc.margin <<- 2^Nvars
                                 }else
                                 {
                                   recalc.margin <<- 2^Nvars
                                 }
                                 last.mutation<<-as.integer(2^(Nvars/2)*0.01)
                                 seed <<- as.integer(runif(n = 1,min = 1,max = 10000))
                                 p.prior <<- runif(n = Nvars, min = 0.5,max = 0.5)
                               }
                               else
                               {
                                 max.cpu <<- as.integer(search.args.list$max.cpu)
                                 objective <<- as.integer(search.args.list$objective)
                                 parallelize <<- search.args.list$parallelize
                                 parallelize.global <<- search.args.list$parallelize.global
                                 parallelize.hyper  <<- search.args.list$parallelize.hyper
                                 p.prior <<-  search.args.list$p.prior
                                 min.N <<- as.integer(search.args.list$min.N)
                                 printable.opt <<- search.args.list$printable.opt
                                 min.N.glob <<- as.integer(search.args.list$min.N.glob)
                                 max.N.glob <<- as.integer(search.args.list$max.N.glob)
                                 switch.type.glob <<- as.integer(search.args.list$switch.type.glob)
                                 min.N.randomize <<- as.integer(search.args.list$min.N.randomize)
                                 max.N.randomize <<- as.integer(search.args.list$max.N.randomize)
                                 type.randomize <<- as.integer(search.args.list$type.randomize)
                                 max.cpu.glob <<- as.integer(search.args.list$max.cpu.glob)
                                 locstop.nd <<-search.args.list$locstop.nd
                                 max.cpu.hyper <<- as.integer(search.args.list$max.cpu.hyper)
                                 save.beta <<- search.args.list$save.beta
                                 aa <<- search.args.list$lambda.a
                                 thin_rate <-search.args.list$thin_rate
                                 cc <<- search.args.list$lambda.c
                                 M.nd <<- as.integer(search.args.list$stepsGreedy)
                                 M.mcmc <<- as.integer(search.args.list$stepsLocMCMC)
                                 SA.param <<- search.args.list$SA.params
                                 fobserved <<- search.args.list$fobserved
                                 switch.type <<- as.integer(search.args.list$fswitch.type)
                                 n.size <<- as.integer(search.args.list$n.size)
                                 LocImprove <<-search.args.list$prior.optimizer.freq
                                 max.N <<- as.integer(search.args.list$max.N)
                                 fparam <<- search.args.list$fparam
                                 fparam.pool <<- search.args.list$fparam
                                 isobsbinary <<- as.array(0:(length(fparam)-1))
                                 p.add <<- as.array(search.args.list$p.add)
                                 recalc.margin <<- search.args.list$recalc.margin
                                 Nvars  <<- as.integer(length(fparam))
                                 seed <<-  search.args.list$seed
                                 max.tree.size<<- as.integer(search.args.list$max.tree.size)
                                 Nvars.max <<- as.integer(search.args.list$Nvars.max)
                                 Nvars.init <<- as.integer(search.args.list$Nvars)
                                 allow_offsprings <<- as.integer(search.args.list$allow_offsprings)
                                 mutation_rate <<- as.integer(search.args.list$mutation_rate)
                                 p.allow.tree <<- search.args.list$p.allow.tree
                                 p.allow.replace <<- search.args.list$p.allow.replace
                                 last.mutation<<-as.integer(search.args.list$last.mutation)
                                 p.nor <<- search.args.list$p.nor
                                 p.and <<- search.args.list$p.and
                                 sigmas<<-search.args.list$sigmas
                                 sigmas.prob<<-search.args.list$sigmas.prob
                                 double.hashing <<- search.args.list$double.hashing
                                 hash.length <<- as.integer(search.args.list$hash.length)
                               }

                             },
                             #transform binary numbers to decimal
                             bittodec.alt = function(bit) #transform a binary vector into a natural number to correspond between vector of solutions and storage array
                             {

                               n<-length(bit)
                               dec <- 0
                               for(i in 1:n)
                               {
                                 j<-n-i
                                 dec <- dec + ((2)^j)*bit[i]
                               }
                               return(dec)
                             },
                             bittodec = function(bit) #transform a binary vector into a natural number to correspond between vector of solutions and storage array
                             {
                               if(!double.hashing){
                                 n<-length(bit)
                                 dec <- 0
                                 for(i in 1:n)
                                 {
                                   j<-n-i
                                   dec <- dec + ((2)^j)*bit[i]
                                 }

                                 return(dec)

                               }
                               else
                               {
                                 if(exists("statistics1"))
                                 {

                                   hash.level<-0
                                   dec<- hashing(bit)+1
                                   #print(dec)
                                   sum.one<- sum(bit)*which.max(na.rm = T,bit) + sum(bit[Nvars -7:Nvars+1])*Nvars
                                   jjj<-1
                                   while(!add.key(dec,bit,sum.one,T))
                                   {
                                     hash.level<-hash.level+1
                                     bit1<- dectobit.alt(2654435761*(dec+97*sum.one+hash.level*36599)+hash.level*59+hash.level)
                                     dec<- hashing(bit1)+1
                                     #jjj<-jjj+1
                                     #print(dec)
                                   }
                                   #print(jjj)
                                   dec <- dec - 1
                                 }
                                 else if(exists("statistics"))
                                 {
                                   hash.level<-0
                                   dec<- hashing(bit)+1
                                   sum.one<- sum(bit)*which.max(na.rm = T,bit) + sum(bit[Nvars -7:Nvars +1])*Nvars
                                   while(!add.key(dec,bit,sum.one,F))
                                   {
                                     hash.level<-hash.level+1
                                     bit1<- dectobit.alt(2654435761*(dec+97*sum.one+hash.level*36599)+hash.level*59+hash.level)
                                     dec<- hashing(bit1)+1
                                     #print(dec)
                                   }
                                   dec <- dec - 1
                                 }
                                 return(dec)
                               }


                             },
                             add.key = function(dec,bit,sum.one,levl)
                             {
                               lb<-length(bit)
                               if(dec>2^hash.length)
                                 return(FALSE)

                               if(levl)
                               {

                                 if(is.na(statistics1[dec,1]))
                                 {
                                   statistics1[dec,16]<-sum.one
                                   hash.keys1[dec,]<-bit
                                   return(TRUE)
                                 }

                                 if(is.na(statistics1[dec,16]))
                                 {
                                   statistics1[dec,16]<-sum.one
                                   hash.keys1[dec,]<-bit#c(array(0,dim = Nvars-lb),bit)
                                   return(TRUE)
                                 }

                                 if(statistics1[dec,16]!=sum.one)
                                   return(FALSE)
                                 i<-1
                                 #print(lb)
                                 lb <- length(bit)
                                 while(i<=lb&&(hash.keys1[dec,Nvars-i+1])==bit[lb - i +1])
                                   i<-i+1
                                 return((i-1)==lb)
                               }
                               else
                               {

                                 if(is.na(statistics[dec,1]))
                                 {
                                   statistics[dec,16]<-sum.one
                                   hash.keys[dec,]<-bit
                                   return(TRUE)
                                 }

                                 if(is.na(statistics[dec,16]))
                                 {
                                   statistics[dec,16]<-sum.one
                                   hash.keys[dec,]<-bit#c(array(0,dim = Nvars-lb),bit)
                                   return(TRUE)
                                 }

                                 if(statistics[dec,16]!=sum.one)
                                   return(FALSE)
                                 i<-1
                                 lb <- length(bit)
                                 while(i<=lb&&(hash.keys[dec,Nvars-i+1])==bit[lb - i +1])
                                   i<-i+1
                                 return((i-1)==lb)
                               }
                             },
                             hashing = function(bit)# a hash function to find where to place the key in the hash
                             {
                               n<-length(bit)
                               if(n<hash.length)
                                 return(bittodec.alt(bit))
                               return(bittodec.alt(bit[(n-hash.length+1):n]))
                             },
                             binlog = function (x) # bitwise logorithm (base 2) computations based on Al Kashi's algorithm
                             {

                               lx <- length(x)
                               tol = -lx
                               y <- 0 # initialise output
                               b <- 0.5 # initialise mantissa


                               # arrange the input into a known range
                               if(lx == 1){
                                 if(x[1] == 0)
                                   return(-Inf)

                                 if(x[1] == 1)
                                   return(0)
                               }
                               # move one bit to the left end
                               # elsewhise

                               if(x[2]==0 && lx==2)
                                 return(1)

                               powto<-2^(-c(1:(lx-1)))
                               float.x<-sum(x[2:lx]*powto)
                               x <- 1 +  float.x
                               y <- lx -1


                               f<-0
                               fb<--1
                               # move one bit to the right end
                               # now x = 1.5
                               # loop until desired tolerance met

                               while(fb > tol)
                               {

                                 x <- x*x

                                 # update the index
                                 if (x >= 2)
                                 {
                                   x <- x/2
                                   y <- y + b
                                   f <- log(exp(f)+exp(fb))
                                 }
                                 # scale for the next bit
                                 b  <- b/2
                                 fb <- (fb-log(2))
                               }
                               return(list(y = y,z = lx -1, f = f))
                             },
                             #transform decimal numbers to binary
                             dectobit = function(dec) #transform a natural number into a binary vector to correspond between vector of solutions and storage array
                             {
                               if(!double.hashing){
                                 if(dec == 0)
                                   return(0)
                                 q<-dec
                                 bin<-NULL
                                 while(q!=0)
                                 {
                                   r<-q/2
                                   q=floor(r)
                                   bin<-c(as.integer(2*(r-q)),bin)
                                 }
                                 return(bin)
                               }

                               return(dehash(dec+1))


                             },
                             dectobit.alt = function(dec) #transform a natural number into a binary vector to correspond between vector of solutions and storage array
                             {

                               if(dec == 0)
                                 return(0)
                               q<-dec
                               bin<-NULL
                               while(q!=0)
                               {
                                 r<-q/2
                                 q=floor(r)
                                 bin<-c(as.integer(2*(r-q)),bin)
                               }
                               return(bin)



                             },
                             dehash=function(dec)
                             {
                               if(exists("hash.keys1"))
                                 return(hash.keys1[dec,])
                               if(exists("hash.keys"))
                                 return(hash.keys[dec,])
                               return(NULL)
                             },
                             #calculate move probabilities
                             calculate.move.logprobabilities = function(varold, varnew, switch.type,min.N,max.N)
                             {
                               if(switch.type == 1) # random size random N(x)
                               {
                                 warning("This option should not be chosen for randomization unless p.add == 0.5 ", call. = FALSE)
                                 min.N = max.N

                                 log.mod.switch.prob <- log(1/(max.N - min.N +1)) # probability of having that many differences
                                 KK<-sum(abs(varold-varnew))

                                 log.mod.switch.prob <- 0 #always the same probabilities for moves within thenighbourhood for p=0.5
                                 log.mod.switchback.prob <- 0
                               } else if(switch.type == 2) #fixed N(x) inverse operator
                               {
                                 if(min.N!=max.N)
                                 {
                                   warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
                                   min.N = max.N
                                 }
                                 log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                 KK<-max.N
                                 log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                 log.mod.switchback.prob <- log.mod.switch.prob
                               }else if(switch.type == 3)  # random sized inverse N(x)
                               {
                                 log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                 KK<-sum(abs(varold-varnew))
                                 log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                 log.mod.switchback.prob <- log.mod.switch.prob
                               }else if(switch.type == 4)  # fixed N(x) for reverse from type 2 swaps
                               {
                                 if(min.N!=max.N)
                                 {
                                   warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
                                   min.N = max.N
                                 }
                                 log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                 KK<-max.N
                                 log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                 log.mod.switchback.prob <- log.mod.switch.prob
                               }else if(switch.type >  4)
                               {

                                 log.mod.switch.prob <- log(x = 1)
                                 log.mod.switchback.prob <-log(x = 1)
                               }
                               return(list(log.switch.forw.prob = log.mod.switch.prob, log.switch.back.prob = log.mod.switchback.prob))
                             },
                             #build a new model to draw from the old one return selection probabilities
                             buildmodel= function(varcur.old,statid, shift.cpu,max.cpu,switch.type,min.N,max.N,changeble.coord = NULL)
                             {
                               vect<- vector(length = max.cpu,mode = "list")
                               shift<-0
                               for(cpu in 1:(max.cpu))
                               {
                                 set.seed(runif(1,1,10000), kind = NULL, normal.kind = NULL)
                                 if(!is.null(varcur.old)&&length(varcur.old==Nvars))
                                 {
                                   varcur.old[which(is.na(varcur.old))]<-1
                                   varcur<-varcur.old
                                 }else
                                 {
                                   varcur<-rbinom(n = (Nvars),size = 1,0.5)
                                   varcur.old<-varcur
                                 }

                                 changeble <-FALSE
                                 if(!is.null(changeble.coord))
                                 {
                                   changeble <- TRUE
                                 }
                                 if(switch.type == 1) # random size random N(x) # try avoiding when doing mcmc rather than optimization
                                 {
                                   log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                   KK<-floor(runif(n=1, min.N, max.N + 0.999999999))
                                   log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                   log.mod.switchback.prob <- log.mod.switch.prob
                                   change.buf <- array(data = 0,dim = Nvars)
                                   if(changeble){
                                     for(ttt in 1:KK)
                                     {
                                       iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
                                       if(changeble.coord[iid]==1)
                                         next
                                       change.buf[iid] = 1
                                       varcur[iid] <-rbinom(n = 1,size = 1,prob =round(p.add[iid],digits = 8))
                                       log.mod.switch.prob = 0#log.mod.switch.prob + log(dbinom(x = varcur[iid],size = 1,prob = p.add[iid])) # this is one of the pathes to get there in general
                                       log.mod.switchback.prob = 0# log.mod.switchback.prob + log(dbinom(x = 1-varcur[iid],size = 1,prob = p.add[iid])) # but this is one of the ways to get back only
                                     }
                                   }else
                                   {
                                     iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
                                     change.buf[iid] = 1
                                     varcur[iid] <-rbinom(n = KK,size = 1,prob = round(p.add,digits = 8))
                                   }
                                 }else if(switch.type == 2) # fixed sized inverse N(x)
                                 {
                                   if(min.N!=max.N)
                                   {
                                     warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
                                     min.N <- max.N
                                   }
                                   log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                   KK<-max.N
                                   log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                   log.mod.switchback.prob <- log.mod.switch.prob

                                   change.buf <- array(data = 0,dim = Nvars)
                                   if(changeble){
                                     for(ttt in 1:KK)
                                     {
                                       iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
                                       if(change.buf[iid]==1 || changeble.coord[iid]==1)
                                       {
                                         KK<-KK+1
                                         next
                                       }
                                       change.buf[iid] = 1
                                       varcur[iid] = 1-varcur[iid]
                                     }
                                   }else
                                   {
                                     iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
                                     change.buf[iid] = 1
                                     varcur[iid] = 1-varcur[iid]
                                   }

                                 }else if(switch.type == 3)  # random sized inverse N(x)
                                 {
                                   log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                   KK<-floor(runif(n=1, min.N, max.N + 0.999999999))
                                   log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                   log.mod.switchback.prob <- log.mod.switch.prob
                                   change.buf <- array(data = 0,dim = Nvars)
                                   if(changeble){
                                     for(ttt in 1:KK)
                                     {
                                       iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
                                       if(change.buf[iid]==1 || changeble.coord[iid]==1)
                                       {
                                         KK<-KK+1
                                         next
                                       }
                                       change.buf[iid] = 1
                                       varcur[iid] = 1-varcur[iid]
                                     }
                                   }else
                                   {
                                     iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
                                     change.buf[iid] = 1
                                     varcur[iid] = 1-varcur[iid]
                                   }
                                 }else if(switch.type == 4)  # fixed N(x) for reverse from type 2 swaps
                                 {
                                   if(min.N!=max.N)
                                   {
                                     warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
                                     min.N <- max.N
                                   }
                                   log.mod.switch.prob <- log(1/(max.N - min.N +1))
                                   KK<-max.N
                                   log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
                                   log.mod.switchback.prob <- log.mod.switch.prob
                                   ids <- which(changeble.coord == 0)
                                   change.buf <- array(data = 0,dim = Nvars)
                                   if(changeble){
                                     for(ttt in KK)
                                     {
                                       iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
                                       if(change.buf[iid]==1 || changeble.coord[iid]==1)
                                       {
                                         KK<-KK+1
                                         next
                                       }
                                       change.buf[iid] = 1
                                       varcur[iid] = 1-varcur[iid]
                                     }}
                                   else
                                   {
                                     iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
                                     change.buf[iid] = 1
                                     varcur[iid] = 1-varcur[iid]
                                   }
                                 }else if(switch.type == 5)
                                 {
                                   change.buf <- array(data = 0,dim = (Nvars))
                                   log.mod.switch.prob<-0
                                   log.mod.switchback.prob <-0
                                   add<-0
                                   if(cpu+shift>Nvars)
                                   {
                                     shift<-1 - cpu
                                     varcur.old<-rep(0,times = Nvars)
                                   }
                                   if(varcur.old[cpu+shift]!=1)
                                   {
                                     varcur<-varcur.old
                                     varcur[cpu+shift]<-1
                                   }else
                                   {
                                     if(cpu+shift<Nvars)
                                       shift<-shift+1
                                     while(varcur.old[cpu+shift]==1)
                                     {
                                       shift<-shift+1
                                       if(cpu+shift>=Nvars)
                                       {
                                         shift <- (Nvars - cpu)
                                         break
                                       }
                                     }
                                     varcur<-varcur.old
                                     varcur[cpu+shift]<-1
                                   }

                                 }else if(switch.type == 6)
                                 {
                                   change.buf <- array(data = 0,dim = (Nvars))
                                   log.mod.switch.prob<-0
                                   log.mod.switchback.prob <-0
                                   add<-0
                                   if(cpu+shift>Nvars)
                                   {
                                     shift<-1 - cpu
                                     varcur.old<-rep(1,times = Nvars)
                                   }

                                   if(varcur.old[cpu+shift]!=0)
                                   {
                                     varcur<-varcur.old
                                     varcur[cpu+shift]<-0
                                   }else
                                   {

                                     if(cpu+shift<Nvars)
                                       shift<-shift+1
                                     while(varcur.old[cpu+shift]==0)
                                     {
                                       shift<-shift+1
                                       if(cpu+shift>=Nvars)
                                       {
                                         shift <- (Nvars - cpu)
                                         break
                                       }
                                     }
                                     varcur<-varcur.old
                                     varcur[cpu+shift]<-0
                                   }

                                 }else if(switch.type == 7)
                                 {
                                   change.buf <- array(data = 0,dim = (Nvars))
                                   log.mod.switch.prob<-0
                                   log.mod.switchback.prob <-0
                                   vec<-dectobit(cpu + shift.cpu)
                                   varcur<-c(array(0,dim = (Nvars -length(vec))),vec) #issues here
                                 }else if(switch.type == 8)
                                 {

                                   log.mod.switch.prob <-0
                                   log.mod.switchback.prob <-0
                                   varcur<-varcur.old
                                   change.buf <- array(data = 0,dim = Nvars)
                                   #if(printable.opt)print("type 8 invoked")
                                   #if(printable.opt)print(varcur)
                                 }else
                                 {

                                   log.mod.switch.prob <- 0
                                   log.mod.switchback.prob <-0
                                   change.buf <- array(data = 1,dim = Nvars)
                                   varcur <-rbinom(n = Nvars,size = 1,prob = round(p.add,digits = 8))
                                   #if(printable.opt)print("type 8 invoked")
                                   #if(printable.opt)print(varcur)
                                 }

                                 #     for(g in 1:max(na.rm = T,isobsbinary))
                                 #     {
                                 #       if(length(varcur[which(isobsbinary == g && varcur %in% c(1,3))])==length(which(isobsbinary == g)))
                                 #         varcur[which(isobsbinary == g && varcur %in% c(1,3))[1]]=0
                                 #       if(length(varcur[which(isobsbinary == g && varcur %in% c(2,3))])==length(which(isobsbinary == g)))
                                 #         varcur[which(isobsbinary == g && varcur %in% c(2,3))[1]]=0
                                 #     }



                                 covobs <- if(fparam[1]=="Const")fparam[which(varcur[-1] == 1)+1]else fparam[which(varcur==1)]


                                 #obsconst<-2*as.integer((varcur[1] %in% c(1,3)) || length(covobs) == 0 || (length(covobs)==length(which(isobsbinary != 0))) ) -1

                                 obsconst<-ifelse(fparam[1]=="Const",2*as.integer((varcur[1])) -1,1)


                                 id<-bittodec(varcur)
                                 id<-id+1



                                 if(ifelse(exists("statistics1"),is.na(statistics1[id,1]),ifelse(exists("statistics"),is.na(statistics[id,1,]),ifelse(exists("hashStat"),!has.key(hash = hashStat,key = paste(varcur,collapse = "")),TRUE))))#||TRUE)
                                 {
                                   # formula <- NULL
                                   # formula <- as.formula(ifelse(length(covobs)>0,(stri_join(stri_flatten(fobserved[1]), " ~ ",obsconst,"+", stri_flatten(covobs, collapse=" + "), latent.formula)),(stri_join(stri_flatten(fobserved[1]), " ~ ",obsconst, latent.formula))))
                                   #
                                   # if(is.null(formula)){
                                   #   formula <- as.formula(ifelse(length(covobs)>0,(stri_join(stri_flatten(fobserved[1]), " ~ ",obsconst,"+", stri_flatten(covobs, collapse=" + "))),(stri_join(stri_flatten(fobserved[1]), " ~ ",obsconst))))
                                   #
                                   # }
                                   formula <- NULL
                                   capture.output({withRestarts(tryCatch(capture.output({formula <- as.formula(paste(paste(fobserved[1]), " ~ ",obsconst,ifelse(length(covobs)>0," + ",""), paste(covobs, collapse=" + "), latent.formula)) })), abort = function(){onerr<-TRUE;fm<-NULL})}) ## not considered currently in RJMCMC, is only valid for model selection
                                   if(is.null(formula)){
                                     formula <- as.formula(paste(paste(fobserved[1]), " ~ ",obsconst,ifelse(length(covobs)>0," + ",""), paste(covobs, collapse=" + ")))

                                   }

                                 }else
                                 {
                                   formula <- NULL
                                 }

                                 vect[[cpu]]<-list(formula = formula, varcur = varcur, statid = statid, changed = change.buf, log.mod.switch.prob = log.mod.switch.prob , log.mod.switchback.prob = log.mod.switchback.prob )

                               }
                               return(vect)
                             },
                             #fit the given model by means of the specified estimator
                             fitmodel = function(model)
                             {
                               if(!is.null(model))
                               {

                                 fm<-NULL
                                 id<-bittodec(model$varcur)

                                 if(is.null(id))
                                   id = 0
                                 id<-id+1

                                 if(exists("statistics1")){
                                   if(is.na(statistics1[id,1]))
                                   {

                                     #if(printable.opt)print("Invoked from EMJMCMC environment")
                                     onerr<-FALSE
                                     #if(printable.opt)print("INLA internal error")
                                     statistics1[id,c(2,3)]<-100000
                                     statistics1[id,1]<--100000
                                     statistics1[id,4:14]<-0

                                     #prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))

                                     #       capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
                                     #       })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
                                     #
                                     capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
                                     })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements



                                     if(!is.null(fm))
                                     {
                                       statistics1[id,2]<-fm$waic[[1]]
                                       statistics1[id,1]<-fm$mlik[[1]]
                                       statistics1[id,3]<-fm$dic[[1]]
                                       if(save.beta)
                                       {


                                         if(fparam[1]=="Const")
                                         {
                                           inxx<-which(model$varcur==1)
                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             statistics1[id,15+inxx]<-fm$summary.fixed$mean
                                         }else
                                         {
                                           inxx<-c(0,which(model$varcur==1))
                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             statistics1[id,16+inxx]<-fm$summary.fixed$mean
                                         }


                                       }

                                       if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
                                       {
                                         g.results[2,1]<-fm$waic[[1]]
                                         g.results[2,2]<-as.integer(id)
                                       }
                                       if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
                                       {
                                         g.results[1,1]<-fm$mlik[[1]]
                                         g.results[1,2]<-as.integer(id)
                                       }

                                       if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
                                       {
                                         g.results[3,1]<-fm$dic[[1]]
                                         g.results[3,2]<-as.integer(id)
                                       }

                                       g.results[4,2] <- g.results[4,2]+1
                                       if(g.results[4,2]%%recalc.margin == 0)
                                       {
                                         p.add <<- as.array(post_proceed_results(statistics1)$p.post)
                                       }

                                     }


                                   }
                                   if(model$statid!=-1)
                                     statistics1[id,model$statid+1]<-statistics1[id,model$statid+1] + 1
                                   g.results[4,1] <- g.results[4,1]+1
                                   return(list(mlik=statistics1[id,1],waic=statistics1[id,2],dic=statistics1[id,3]))
                                 }else  if(exists("statistics")){
                                   if(is.na(statistics[id,1]))
                                   {
                                     #if(printable.opt)print("Invoked from EMJMCMC SUB environment")
                                     onerr<-FALSE
                                     #if(printable.opt)print("INLA internal error")
                                     statistics[id,c(2,3)]<-100000
                                     statistics[id,1]<--100000
                                     statistics[id,4:14]<-0
                                     #prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))

                                     #       capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
                                     #       })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
                                     #
                                     capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
                                     })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements



                                     if(!is.null(fm))
                                     {
                                       statistics[id,2]<-fm$waic[[1]]
                                       statistics[id,1]<-fm$mlik[[1]]
                                       statistics[id,3]<-fm$dic[[1]]
                                       if(save.beta)
                                       {


                                         if(fparam[1]=="Const")
                                         {
                                           inxx<-which(model$varcur==1)
                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             statistics[id,15+inxx]<-fm$summary.fixed$mean
                                         }else
                                         {
                                           inxx<-c(0,which(model$varcur==1))
                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             statistics[id,16+inxx]<-fm$summary.fixed$mean
                                         }


                                       }

                                       if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
                                       {
                                         g.results[2,1]<-fm$waic[[1]]
                                         g.results[2,2]<-as.integer(id)
                                       }
                                       if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
                                       {
                                         g.results[1,1]<-fm$mlik[[1]]
                                         g.results[1,2]<-as.integer(id)
                                       }

                                       if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
                                       {
                                         g.results[3,1]<-fm$dic[[1]]
                                         g.results[3,2]<-as.integer(id)
                                       }

                                       g.results[4,2] <- g.results[4,2]+1
                                       if(g.results[4,2]%%recalc.margin == 0)
                                       {
                                         proceeeded <- post_proceed_results(statistics)
                                         p.add <<- as.array(proceeeded$p.post)
                                         #g.results[4,2] <-
                                       }
                                     }


                                   }
                                   if(model$statid!=-1)
                                     statistics[id,model$statid+1]<-statistics[id,model$statid+1] + 1
                                   g.results[4,1] <- g.results[4,1]+1
                                   return(list(mlik=statistics[id,1],waic=statistics[id,2],dic=statistics[id,3]))

                                 }else if(exists("hashStat")){
                                   if(Nvars.max>Nvars)
                                     idd<- as.character(paste(c(model$varcur,array(0,Nvars.max-Nvars)),collapse = ""))
                                   else
                                     idd<- as.character(paste(c(model$varcur),collapse = ""))

                                   if(!has.key(key = idd,hash = hashStat))
                                   {
                                     #if(printable.opt)print("Invoked from EMJMCMC hash table environment")
                                     onerr<-FALSE
                                     #if(printable.opt)print("INLA internal error")

                                     #prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))

                                     #       capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
                                     #       })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
                                     #
                                     capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
                                     })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements

                                     if(!save.beta)
                                     {
                                       hashBuf<-array(data = NA,dim = 3)
                                     }else
                                     {
                                       if(allow_offsprings==0||Nvars.max<Nvars)
                                       {
                                         if(fparam[1]=="Const")
                                         {
                                           linx<-Nvars
                                           inxx<-which(model$varcur==1)
                                         }else
                                         {
                                           linx<-Nvars + 1
                                           inxx<-c(0,which(model$varcur==1))
                                         }
                                       }else
                                       {
                                         if(fparam[1]=="Const")
                                         {
                                           linx<-Nvars.max
                                           inxx<-which(model$varcur==1)
                                         }else
                                         {
                                           linx<-Nvars.max + 1
                                           inxx<-c(0,which(model$varcur==1))
                                         }
                                       }

                                       hashBuf<-array(data = NA,dim = 3 + linx)
                                     }
                                     hashBuf[c(2,3)]<-100000
                                     hashBuf[1]<- -100000


                                     if(!is.null(fm))
                                     {
                                       hashBuf[2]<-fm$waic[[1]]
                                       hashBuf[1]<-fm$mlik[[1]]
                                       hashBuf[3]<-fm$dic[[1]]

                                       if(save.beta)
                                       {


                                         if(fparam[1]=="Const")
                                         {

                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             hashBuf[3+inxx]<-fm$summary.fixed$mean
                                         }else
                                         {

                                           if(length(inxx)==length(fm$summary.fixed$mean))
                                             hashBuf[4+inxx]<-fm$summary.fixed$mean
                                         }


                                       }

                                       hashStat[idd] <- hashBuf
                                       #                                                          if(id>1)
                                       #                                                          {
                                       #                                                            inxx<-which(model$varcur==1)
                                       #                                                            if(length(inxx)==length(fm$summary.fixed$mean))
                                       #                                                              statistics[id,14+inxx]<-fm$summary.fixed$mean
                                       #                                                          }
                                       if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
                                       {
                                         g.results[2,1]<-fm$waic[[1]]
                                         g.results[2,2]<-(id)
                                       }
                                       if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
                                       {
                                         g.results[1,1]<-fm$mlik[[1]]
                                         g.results[1,2]<-(id)
                                       }

                                       if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
                                       {
                                         g.results[3,1]<-fm$dic[[1]]
                                         g.results[3,2]<-(id)
                                       }

                                       g.results[4,2] <- g.results[4,2]+1

                                     }


                                   }
                                   g.results[4,1] <- g.results[4,1]+1
                                   if(has.key(hash = hashStat,key=idd))
                                     hasRes<- values(hashStat[idd])
                                   else
                                     hasRes<-c(-10000,10000,10000)
                                   if(g.results[4,2]%%recalc.margin == 0)
                                   {
                                     proceeeded <- post_proceed_results_hash(hashStat)
                                     p.add <<- as.array(proceeeded$p.post)
                                     #g.results[4,2] <-
                                   }
                                   return(list(mlik=hasRes[1],waic=hasRes[2],dic=hasRes[3]))

                                 }else
                                 {
                                   capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
                                   })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements

                                   if(!is.null(fm)){

                                     if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
                                     {
                                       g.results[2,1]<-fm$waic[[1]]
                                       g.results[2,2]<-(id)
                                     }
                                     if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
                                     {
                                       g.results[1,1]<-fm$mlik[[1]]
                                       g.results[1,2]<-(id)
                                     }

                                     if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
                                     {
                                       g.results[3,1]<-fm$dic[[1]]
                                       g.results[3,2]<-(id)
                                     }
                                     g.results[4,1] <- g.results[4,1]+1
                                     g.results[4,2] <- g.results[4,2]+1
                                     return(list(mlik=fm$mlik[[1]],waic=fm$waic[[1]],dic=fm$dic[[1]]))
                                   }
                                   else
                                   {
                                     g.results[4,1] <- g.results[4,1]+1
                                     return(list(mlik=-Inf,waic=Inf,dic=Inf))
                                   }
                                 }
                               }
                               g.results[4,1] <- g.results[4,1]+1
                               g.results[4,2] <- g.results[4,2]+1
                               return(list(mlik=-Inf,waic=Inf,dic=Inf))
                             },
                             #lambda function for mtmcmc
                             lambda = function(c,alpha,g1,g2,g.domain.pos) # simmetric choice driving function
                             {
                               if((c!=0))
                               {
                                 res<-((((1/c)*(1+(g1+g2)*(g.domain.pos))/(1-(g1+g2)*(1-g.domain.pos)))^alpha))
                                 #!#if(printable.opt)print(res)
                                 return(res)
                               }else
                               {
                                 #!#if(printable.opt)print(1)
                                 return(1)
                               }
                             },
                             #norm between probabilities vectors
                             normprob=function(p1,p2)
                             {
                               nn<-abs(1-sum((p1+0.1)/(p2+0.1))/length(p1))
                               if(is.na(nn))
                                 nn<-Inf
                               return(nn)
                             },
                             #local mcmc procedure
                             learnlocalMCMC=function(model)
                             {
                               M<-M.mcmc
                               mlikcur<- model$mlikcur
                               waiccur<-model$waiccur
                               varcand<-model$varcur
                               varcur<-model$varcur
                               varcurb<-model$varcur
                               varglob<-model$varcur

                               modglob<-NULL
                               fm<-NULL
                               fmb<-NULL

                               # estimate large jump in a reverse move
                               if(model$reverse || is.infinite(mlikcur))
                               {
                                 vectbg<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
                                 if(!is.null(vectbg[[1]]$formula))
                                 {
                                   bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
                                   waiccur<-bgmod[[1]]$waic
                                   mlikcur<-bgmod[[1]]$mlik
                                 }
                                 else if(exists("statistics1"))
                                 {
                                   iidd<-bittodec(varcur)+1
                                   waiccur<-statistics1[iidd,2]
                                   mlikcur<-statistics1[iidd,1]
                                 }else if(exists("hashStat"))
                                 {
                                   iidd<-paste(varcur,collapse = "")
                                   waiccur<-values(hashStat[iidd])[2]
                                   mlikcur<-values(hashStat[iidd])[1]
                                 }



                               }


                               if(printable.opt)print(paste("Begin with ", mlikcur))

                               mlikglob<- mlikcur
                               mlikcand<- mlikcur
                               waiccand<- waiccur
                               waicglob<- waiccur
                               waiccur<-  waiccur

                               for(m in 1:M)
                               {
                                 withRestarts(tryCatch({

                                   # statistics <- describe(statistics)
                                   vect<-buildmodel(varcur.old = varcur,statid = model$statid,max.cpu = max.cpu,switch.type=switch.type, min.N = min.N, max.N = max.N)
                                   cluster<-TRUE
                                   flag1<-0

                                   for(mod_id in 1:max.cpu)
                                   {
                                     if(is.null(vect[[mod_id]]$formula))
                                     {
                                       flag1<-flag1+1
                                     }

                                   }

                                   #flag1<-sum(is.null(vect[[]]$formula))

                                   if(flag1==max.cpu)
                                   {
                                     cluster<-FALSE
                                     if(printable.opt)print("!!!!MTMCMC models already estimated!!!!")
                                   }else
                                   {

                                     res.par <- parallelize(X = vect,FUN = .self$fitmodel)

                                   }
                                   p.select.y <- array(data = 0, dim = max.cpu)
                                   for(mod_id in 1:max.cpu)
                                   {
                                     if(cluster)
                                     {
                                       fm<-res.par[[mod_id]]

                                       if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
                                       {
                                         varcand<-varcurb
                                         if(printable.opt)print("locMTMCMC Model Fit Error!?")
                                         next
                                       }
                                     }

                                     varcand<-vect[[mod_id]]$varcur

                                     if(cluster)
                                     {
                                       waiccand<-res.par[[mod_id]]$waic
                                       mlikcand<-res.par[[mod_id]]$mlik
                                     }else if(exists("statistics1"))
                                     {
                                       iidd<-bittodec(varcand)+1
                                       waiccand<-statistics1[iidd,2]
                                       mlikcand<-statistics1[iidd,1]
                                     }else if(exists("hashStat"))
                                     {
                                       iidd<-paste(varcand,collapse = "")
                                       waiccand<-values(hashStat[iidd])[2]
                                       mlikcand<-values(hashStat[iidd])[1]
                                     }

                                     if((mlikcand>mlikglob)) #update the parameter of interest
                                     {
                                       if(printable.opt)print(paste("locMTMCMC update waic.glob = ", waiccand))
                                       if(printable.opt)print(paste("locMTMCMC update waic.glob.mlik = ",  mlikglob))
                                       mlikglob<-mlikcand
                                       waicglob<-waiccand
                                       varglob<-varcand
                                       if(cluster)
                                         modglob<-fm
                                     }


                                     g1 <- waiccur

                                     if(waiccur == Inf)
                                     {
                                       g1 = 1
                                     }

                                     p.select.y[mod_id]<-(mlikcand + vect[[mod_id]]$log.mod.switchback.prob+log(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos =  FALSE))) # correct for different criteria later

                                     if(is.na(p.select.y[mod_id]))
                                       p.select.y[mod_id] <- 0
                                     if(is.infinite(p.select.y[mod_id]) || p.select.y[mod_id]>100000000)
                                     {
                                       #if(printable.opt)print(paste("very large log.w.y detected ",p.select.y[mod_id]))
                                       p.select.y[mod_id] <- 100000000
                                     }

                                   }

                                   max.p.select.y <- max(na.rm = T,p.select.y)
                                   p.select.y<-p.select.y-max.p.select.y

                                   #if(printable.opt)print(paste("max log.w.y is ",max.p.select.y,"normilized log.w.n.y is ", paste(p.select.y,collapse = ", ")))


                                   ID<-sample(x = max.cpu,size = 1,prob = exp(p.select.y))

                                   if(printable.opt)print(paste("cand ",ID," selected"))

                                   varcand<-vect[[ID]]$varcur

                                   if(cluster)
                                   {
                                     waiccand<-res.par[[ID]]$waic
                                     mlikcand<-res.par[[ID]]$mlik
                                   }else if(exists("statistics1"))
                                   {
                                     iidd<-bittodec(varcand)+1
                                     waiccand<-statistics1[iidd,2]
                                     mlikcand<-statistics1[iidd,1]
                                   }else if(exists("hashStat"))
                                   {
                                     iidd<-paste(varcand,collapse = "")
                                     waiccand<-values(hashStat[iidd])[2]
                                     mlikcand<-values(hashStat[iidd])[1]
                                   }

                                   #p.Q.cand<- p.select.y[ID]/sum(p.select.y)

                                   if(printable.opt)print("do reverse step")

                                   p.select.z <- array(data = 0.01, dim = max.cpu)


                                   if(max.cpu!=1)
                                   {
                                     vect1<-buildmodel(max.cpu = max.cpu -1,varcur.old = varcand,statid = model$statid,switch.type = switch.type, min.N = min.N, max.N = max.N)

                                     cluster<-TRUE

                                     flag1<-0

                                     for(mod_id in 1:(max.cpu-1))
                                     {
                                       if(is.null(vect1[[mod_id]]$formula))
                                       {
                                         flag1<-flag1+1
                                       }

                                     }

                                     if(flag1==(max.cpu-1))
                                     {
                                       cluster<-FALSE
                                       if(printable.opt)print("!!!!MTMCMC reverse models already estimated!!!!")
                                     }else
                                     {
                                       res.par.back <- parallelize(X = vect1,FUN = .self$fitmodel)
                                     }

                                     for(mod_id in 1:(max.cpu-1))
                                     {

                                       if(cluster)
                                       {
                                         if(is.null(fm)&&(is.na(res.par.back[[mod_id]]$waic)))
                                         {
                                           if(printable.opt)print("locMTMCMC Model Fit Error!?")
                                           next
                                         }
                                       }

                                       varcand.b<-vect1[[mod_id]]$varcur

                                       if(cluster)
                                       {
                                         waiccand.b<-res.par.back[[mod_id]]$waic
                                         mlikcand.b<-res.par.back[[mod_id]]$mlik
                                       }else if(exists("statistics1"))
                                       {
                                         iidd<-bittodec(varcand.b)+1
                                         waiccand.b<-statistics1[iidd,2]
                                         mlikcand.b<-statistics1[iidd,1]
                                       }else if(exists("hashStat"))
                                       {
                                         iidd<-paste(varcand.b,collapse = "")
                                         waiccand.b<-values(hashStat[iidd])[2]
                                         mlikcand.b<-values(hashStat[iidd])[1]
                                       }

                                       if((mlikcand.b>mlikglob))
                                       {
                                         if(printable.opt)print(paste("locMTMCMC update waic.glob = ", waiccand.b))
                                         if(printable.opt)print(paste("locMTMCMC update waic.glob.mlik = ", mlikcand.b))
                                         mlikglob<-mlikcand.b
                                         waicglob<-waiccand.b
                                         varglob<-varcand.b
                                         if(cluster)
                                           modglob<-fm
                                       }

                                       g1 = waiccand

                                       if(waiccand == Inf)
                                       {
                                         g1 = 1
                                       }

                                       p.select.z[mod_id]<-(mlikcand.b+vect1[[mod_id]]$log.mod.switchback.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand.b,g.domain.pos = FALSE))) # correct for different criteria later

                                       if(is.na(p.select.z[mod_id]))
                                         p.select.z[mod_id]=0
                                       if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
                                       {
                                         #if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
                                         p.select.z[mod_id] <- 100000000
                                       }
                                     }
                                   }

                                   if( waiccur == Inf)
                                   {
                                     g1 = 1
                                   }
                                   p.select.z[max.cpu] <- (mlikcur+vect[[ID]]$log.mod.switch.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos = FALSE)))

                                   if(is.na(p.select.z[mod_id]))
                                     p.select.z[mod_id]=0
                                   if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
                                   {
                                     #if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
                                     p.select.z[mod_id] <- 100000000
                                   }

                                   max.p.select.z <- max(na.rm = T,p.select.z)
                                   p.select.z<-p.select.z-max.p.select.z

                                   if(printable.opt)print(paste("max log.w.z is ",max.p.select.z,"normilized log.w.n.z is ", paste(p.select.z,collapse = ", ")))

                                   if(log(runif(n = 1,min = 0,max = 1)) < sum(na.rm = T,sum(na.rm = T,log(sum(na.rm = T,exp(p.select.y))),-log(sum(na.rm = T,exp(p.select.z))), max.p.select.y, - max.p.select.z )))
                                   {
                                     mlikcur<-mlikcand
                                     if(printable.opt)print(paste("locMTMCMC update ratcur = ", mlikcand))
                                     if(printable.opt)print(paste("locMTMCMC accept move with ", waiccand))
                                     varcur<-varcand
                                     waiccur<-waiccand

                                   }

                                 }),abort = function(){fm<-fmb;closeAllConnections();options(error=traceback);  onerr<-TRUE})
                               }

                               #if(printable.opt)print("FINISH LOCAL MTMCMC")

                               #!#if(printable.opt)print(points)

                               if(model$reverse == FALSE)
                               {
                                 if(is.null(varcur))
                                 {
                                   #if(printable.opt)print("No moves acceoted in the procedure")
                                   varcur<-varcand
                                   waiccur<-waiccand
                                   varcur<-varcand
                                 }

                                 vect<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)


                                 varcur<-vect[[1]]$varcur
                                 #if(printable.opt)print(varcur)

                                 cluster<-TRUE


                                 if(is.null(vect[[1]]$formula))
                                 {
                                   cluster<-FALSE
                                   if(printable.opt)print("!!!!MTMCMC reverse models already estimated!!!!")
                                 }else
                                 {
                                   mod<-fitmodel(vect[[1]])
                                 }

                                 if(cluster)
                                 {
                                   waiccur<-mod$waic
                                   mlikcur<-mod$mlik
                                 }else if(exists("statistics1"))
                                 {
                                   iidd<-bittodec(varcur)+1
                                   waiccur<-statistics1[iidd,2]
                                   mlikcur<-statistics1[iidd,1]
                                 }else if(exists("hashStat"))
                                 {
                                   iidd<-paste(varcur,collapse = "")
                                   waiccur<-values(hashStat[iidd])[2]
                                   mlikcur<-values(hashStat[iidd])[1]
                                 }
                                 # incorporate what happens for the backward optimization

                                 model.prob<-vect[[1]]$log.mod.switch.prob
                                 model.prob.fix<-vect[[1]]$log.mod.switchback.prob

                               }else  # incorporate what happens for the reverse move
                               {
                                 if(is.null(varcur))
                                 {
                                   #if(printable.opt)print("No moves acceoted in the reverse procedure")
                                   varcur<-varcand
                                   waiccur<-waiccand
                                 }
                                 model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = model$varold,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)
                                 model.prob<-model.probs$log.switch.forw.prob
                                 model.prob.fix<-model.probs$log.switch.back.prob
                               }

                               if(is.null(varcur))
                               {
                                 #if(printable.opt)print("NO VARCUR OBTAINED")
                                 varcur<-i
                                 waiccur<-Inf
                                 varcur<-model$varold
                               }

                               return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
                             },
                             #local simulated annealing optimization
                             learnlocalSA=function(model)
                             {
                               t.min<-SA.param$t.min
                               t<-SA.param$t.init
                               dt<-SA.param$dt
                               M<-SA.param$M
                               varcand<-model$varcur
                               varcurb<-model$varcur
                               varcur<-model$varcur
                               varglob<-model$varcur
                               varglob<-NULL
                               modglob<-NULL
                               probcur<-1
                               probrevcur<-1
                               first.prob <- 1
                               fm<-NULL
                               fmb<-NULL

                               mlikcur<- model$mlikcur
                               waiccur<-model$waiccur
                               # estimate large jump in a reverse move
                               # estimate large jump in a reverse move
                               if((model$reverse && !model$sa2) || is.infinite(mlikcur))
                               {
                                 vectbg<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
                                 if(!is.null(vectbg[[1]]$formula))
                                 {
                                   bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
                                   waiccur<-bgmod[[1]]$waic
                                   mlikcur<-bgmod[[1]]$mlik
                                 }
                                 else if(exists("statistics1"))
                                 {
                                   iidd<-bittodec(varcur)+1
                                   waiccur<-statistics1[iidd,2]
                                   mlikcur<-statistics1[iidd,1]
                                 }else if(exists("hashStat"))
                                 {
                                   iidd<-paste(varcur,collapse = "")
                                   waiccur<-values(hashStat[iidd])[2]
                                   mlikcur<-values(hashStat[iidd])[1]
                                 }

                               }

                               if(printable.opt)print(paste("Begin with ", mlikcur))
                               mlikglob<- mlikcur
                               mlikcand<- mlikcur
                               waiccand<- waiccur
                               waicglob<- waiccur
                               waiccur<-  waiccur

                               while(t>t.min)
                               {
                                 if(printable.opt)print(paste("anneal to ",t))
                                 t.new<-t*exp(-dt)
                                 if(model$reverse == TRUE && t.new <= t.min)
                                 {
                                   M<-M -1
                                 }
                                 for(m in 1:M)
                                 {
                                   withRestarts(tryCatch({

                                     mmax.cpu = max.cpu
                                     if(model$switch.type == 5)
                                     {
                                       mmax.cpu = Nvars - sum(varcur)
                                     }else if(model$switch.type == 6)
                                     {
                                       mmax.cpu = sum(varcur)
                                     }
                                     if(mmax.cpu == 0)
                                       mmax.cpu = 1

                                     vect<-buildmodel(max.cpu = mmax.cpu,varcur.old = varcur,statid = model$statid,switch.type=model$switch.type, min.N = min.N,max.N = max.N)
                                     cluster<-TRUE
                                     flag1<-0
                                     for(mod_id in 1:mmax.cpu)
                                     {
                                       if(is.null(vect[[mod_id]]$formula))
                                       {
                                         flag1<-flag1+1
                                       }

                                     }
                                     if(flag1==mmax.cpu)
                                     {
                                       cluster<-FALSE
                                       if(printable.opt)print("!!!!SA Models already estimated!!!!")
                                     }else
                                     {
                                       res.par <- parallelize(X = vect,FUN = .self$fitmodel)
                                     }
                                     for(mod_id in 1:mmax.cpu)
                                     {
                                       if(cluster)
                                       {
                                         fmb<-fm
                                         fm<-res.par[[mod_id]]
                                         waiccand<-Inf
                                         if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
                                         {
                                           varcand<-varcurb
                                           if(printable.opt)print("SA Model Fit Error!?")
                                           next
                                         }
                                       }

                                       varcand<-vect[[mod_id]]$varcur
                                       if(cluster)
                                       {
                                         waiccand<-res.par[[mod_id]]$waic
                                         mlikcand<-res.par[[mod_id]]$mlik
                                       }else if(exists("statistics1"))
                                       {
                                         iidd<-bittodec(varcand)+1
                                         waiccand<-statistics1[iidd,2]
                                         mlikcand<-statistics1[iidd,1]
                                       }else if(exists("hashStat"))
                                       {
                                         iidd<-paste(varcand,collapse = "")
                                         waiccand<-values(hashStat[iidd])[2]
                                         mlikcand<-values(hashStat[iidd])[1]
                                       }

                                       if(objective==0)
                                       {
                                         objcand<-waiccand
                                         objcur<-waiccur
                                         objglob<-waicglob
                                       }else
                                       {
                                         objcand<- -mlikcand
                                         objcur<-  -mlikcur
                                         objglob<- -mlikglob
                                       }

                                       if(t == SA.param$t.init && mod_id == 1 && m == 2)
                                       {
                                         delta<-objcand - objcur
                                         first.prob <- vect[[mod_id]]$log.mod.switchback.prob + log(punif(q = exp(x = delta/t),min = 0,max = 1))
                                       }

                                       if(objcand<objcur)
                                       {
                                         if(printable.opt)print(paste("SA accept move with ", objcand))
                                         waiccur<-waiccand
                                         varcur<-varcand
                                         mlikcur<-mlikcand
                                         if(objcand<objglob)
                                         {
                                           waicglob<-waiccand
                                           varglob<-varcand
                                           mlikglob<-mlikcand
                                           if(cluster)
                                             modglob<-fm

                                           if(printable.opt)print(paste("SA update global optima with", objcand))
                                         }
                                       }else
                                       {
                                         delta<-objcand - objcur
                                         if(runif(n = 1,min = 0,max = 1) <= sum(na.rm = T, exp(x = -delta/t)))
                                         {

                                           model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = varcand,switch.type = model$switch.type,min.N = min.N,max.N = max.N)
                                           probcur<-model.probs$log.switch.forw.prob
                                           probrevcur<-model.probs$log.switch.back.prob
                                           waiccur<-waiccand
                                           varcur<-varcand
                                           mlikcur<-mlikcand
                                           if(printable.opt)print(paste("SA accept move with ", objcand))
                                         }
                                       }

                                     }

                                   }),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback);  onerr<-TRUE})
                                 }
                                 t<-t.new
                               }
                               t<-t/exp(-dt)

                               if(model$reverse == FALSE)
                               {
                                 model.prob<-log(punif(q = sum(na.rm = T, exp(x = -delta/t)),min = 0,max = 1)) +  probcur # log(P(Mk,Mk-1))
                                 model.prob.fix<-log(punif(q = exp(x = delta/t),min = 0,max = 1)) + probrevcur # log(P(Mk-1,Mk))

                                 if(model$sa2 == TRUE)
                                 {
                                   model.prob.fix<-model.prob.fix + first.prob # correcting for the term for local improvements of type 3.
                                 }
                               }else  # incorporate what happens for the reverse move
                               {

                                 if(is.null(varcur))
                                 {
                                   if(printable.opt)print("No moves accepted in the reverse procedure")
                                   varcur<-model$varcur
                                   objcur<-model$objcur
                                 }

                                 delta<-objcur-model$objold


                                 model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = model$varold, switch.type = model$switch.type,min.N = min.N,max.N = max.N)
                                 model.prob<-punif(q = sum(na.rm = T, exp(x = -delta/t)),min = 0,max = 1,log.p = TRUE) +  model.probs$log.switch.forw.prob
                                 model.prob.fix<-punif(q = exp(x = delta/t),min = 0,max = 1,log.p = TRUE) + model.probs$log.switch.back.prob

                                 if(model.prob==-Inf)
                                 {
                                   model.prob<- -100000000
                                 }
                                 if(model.prob.fix==-Inf)
                                 {
                                   model.prob.fix<- -100000000
                                 }
                               }

                               return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
                             },
                             #forward selection procedure
                             forward_selection=function(model)
                             {
                               if(printable.opt)print("begin forward selection procedure")
                               varcand<-model$varcur
                               varcurb<-model$varcur
                               varglob<-varcand
                               mlikglob<- model$mlikcur
                               mlikcur<- model$mlikcur
                               waiccand<- model$waiccur
                               waicglob<-model$waiccur
                               waiccur<-model$waiccur
                               waiccurb<-model$waiccur
                               varglob<-NULL
                               modglob<-NULL


                               fm<-NULL
                               fmb<-NULL

                               ub<-bittodec(array(1,length(varcurb)))
                               layer<-length(which(varcurb == 0))

                               mlikcur<- model$mlikcur
                               waiccur<-model$waiccur
                               # estimate large jump in a reverse move
                               # estimate large jump in a reverse move
                               if(is.infinite(mlikcur))
                               {
                                 vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
                                 if(!is.null(vectbg[[1]]$formula))
                                 {
                                   bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
                                   waiccur<-bgmod[[1]]$waic
                                   mlikcur<-bgmod[[1]]$mlik
                                 }
                                 else if(exists("statistics1"))
                                 {
                                   iidd<-bittodec(varcand)+1
                                   waiccur<-statistics1[iidd,2]
                                   mlikcur<-statistics1[iidd,1]
                                 }
                                 else if(exists("hashStat"))
                                 {
                                   iidd<-paste(varcand,collapse = "")
                                   waiccur<-values(hashStat[iidd])[2]
                                   mlikcur<-values(hashStat[iidd])[1]
                                 }
                                 if(!is.na(mlikcur) &&  !is.na(waiccur) )
                                 {
                                   mlikglob<- mlikcur
                                   mlikcur<-  mlikcur
                                   waiccand<- waiccur
                                   waicglob<- waiccur
                                   waiccur<-  waiccur
                                   waiccurb<- waiccur
                                 }
                               }


                               while(layer>0)
                               {
                                 withRestarts(tryCatch({

                                   if(printable.opt)print(paste("proceed with layer",layer))
                                   if(printable.opt)print(paste("current solution is",as.character(varcand)))

                                   vect<-buildmodel(max.cpu = layer,varcur.old = varcurb,statid = model$statid, switch.type = 5,min.N = min.N, max.N = max.N)

                                   if(printable.opt)print(paste("finish preparing models at layer",layer))

                                   cluster<-TRUE
                                   flag1<-0
                                   for(mod_id in 1:layer)
                                   {
                                     if(is.null(vect[[mod_id]]$formula))
                                     {
                                       flag1<-flag1+1
                                     }

                                   }
                                   if(flag1==layer)
                                   {
                                     cluster<-FALSE
                                     if(printable.opt)print("!!!!forward Models already estimated!!!!")
                                   }else
                                   {
                                     res.par <- parallelize(X = vect,FUN = .self$fitmodel)
                                   }

                                   if(printable.opt)print(paste("end forward optimizing at layer",layer))
                                   for(mod_id in 1:layer)
                                   {
                                     if(cluster){
                                       fmb<-fm
                                       fm<-res.par[[mod_id]]
                                       waiccand<-Inf
                                       if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
                                       {
                                         varcand<-varcurb
                                         if(printable.opt)print("forward Model Fit Error!?")
                                         next
                                       }
                                     }

                                     varcand<-vect[[mod_id]]$varcur
                                     if(cluster)
                                     {
                                       waiccand<-res.par[[mod_id]]$waic
                                       mlikcand<-res.par[[mod_id]]$mlik
                                     }else if(exists("statistics1"))
                                     {
                                       iidd<-bittodec(varcand)+1
                                       waiccand<-statistics1[iidd,2]
                                       mlikcand<-statistics1[iidd,1]
                                     }else if(exists("hashStat"))
                                     {
                                       iidd<-paste(varcand,collapse = "")
                                       waiccand<-values(hashStat[iidd])[2]
                                       mlikcand<-values(hashStat[iidd])[1]
                                     }

                                     if(objective==0)
                                     {
                                       objcand<-waiccand
                                       objcur<-waiccur
                                       objglob<-waicglob
                                     }else
                                     {
                                       objcand<- -mlikcand
                                       objcur<-  -mlikcur
                                       objglob<- -mlikglob
                                     }


                                     if(objcand<=objcur || mod_id ==1)
                                     {
                                       if(printable.opt)print(paste("forward accept with ", objcand))
                                       objcur<-objcand
                                       varcurb<-varcand
                                       waiccur<-waiccand
                                       varcur<-varcand
                                       mlikcur<-mlikcand

                                       if(objcur<objglob)
                                       {
                                         objglob<-objcur
                                         waicglob<-waiccand
                                         varglob<-varcand
                                         mlikglob<-mlikcand
                                         if(!is.null(fm))
                                           modglob<-fm

                                         if(printable.opt)print(paste("forward global optima with ", objcand))
                                       }
                                     }
                                     #                                                         else
                                     #                                                          {
                                     #                                                            if(waiccand<=waiccur)
                                     #                                                            {
                                     #                                                              waiccur<-waiccand
                                     #                                                              varcur<-varcand
                                     #                                                            }
                                     #
                                     #                                                          }

                                   }
                                   if(objcur!=objglob)
                                   {
                                     if(model$locstop)
                                     {
                                       break
                                     }else
                                     {
                                       varcurb<-varcur
                                     }
                                   }

                                 }),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback);  onerr<-TRUE})


                                 layer<-layer-1
                               }

                               model.prob<-1


                               model.prob.fix<-1


                               return(list(varcur = varglob, waiccur = waicglob, mlikcur = mlikglob, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
                             },
                             #backward selection procedure
                             backward_selection=function(model)
                             {
                               #if(printable.opt)print("begin backward selection procedure")
                               varcand<-model$varcur
                               varcurb<-model$varcur
                               varglob<-varcand
                               mlikglob<- model$mlikcur
                               mlikcur<- model$mlikcur
                               waiccur<- model$waiccur
                               waicglob<-model$waiccur
                               varglob<-NULL
                               modglob<-NULL
                               waiccurb<- model$waiccur

                               fm<-NULL
                               fmb<-NULL

                               ub<-bittodec(array(1,length(varcurb)))
                               layer<-length(which(varcurb == 1))

                               if(is.infinite(mlikcur))
                               {
                                 vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
                                 if(!is.null(vectbg[[1]]$formula))
                                 {
                                   bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
                                   waiccur<-bgmod[[1]]$waic
                                   mlikcur<-bgmod[[1]]$mlik
                                 }
                                 else if(exists("statistics1"))
                                 {
                                   iidd<-bittodec(varcand)+1
                                   waiccur<-statistics1[iidd,2]
                                   mlikcur<-statistics1[iidd,1]
                                 }else if(exists("hashStat"))
                                 {
                                   iidd<-paste(varcand,collapse = "")
                                   waiccur<-values(hashStat[iidd])[2]
                                   mlikcur<-values(hashStat[iidd])[1]
                                 }
                                 if(!is.na(mlikcur) &&  !is.na(waiccur) )
                                 {
                                   mlikglob<- mlikcur
                                   mlikcur<-  mlikcur
                                   waiccand<- waiccur
                                   waicglob<- waiccur
                                   waiccur<-  waiccur
                                   waiccurb<- waiccur
                                 }
                               }


                               while(layer>0)
                               {
                                 withRestarts(tryCatch({

                                   if(printable.opt)print(paste("backward proceed with layer",layer))
                                   if(printable.opt)print(paste("current backward solution is",as.character(varcand)))
                                   vect<-buildmodel(max.cpu = layer,varcur.old = varcurb,statid = model$statid, switch.type = 6,min.N = min.N, max.N = max.N)

                                   if(printable.opt)print(paste("finish backward preparing models at layer",layer))

                                   cluster<-TRUE
                                   flag1<-0
                                   for(mod_id in 1:layer)
                                   {
                                     if(is.null(vect[[mod_id]]$formula))
                                     {
                                       flag1<-flag1+1
                                     }

                                   }
                                   if(flag1==layer)
                                   {
                                     cluster<-FALSE
                                     if(printable.opt)print("!!!!backward Models already estimated!!!!")
                                   }else
                                   {
                                     res.par <- parallelize(X = vect,FUN = .self$fitmodel)
                                   }
                                   if(printable.opt)print(paste("end backward optimizing at layer",layer))

                                   for(mod_id in 1:layer)
                                   {
                                     if(cluster){
                                       fmb<-fm
                                       fm<-res.par[[mod_id]]
                                       waiccand<-Inf
                                       if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
                                       {
                                         varcand<-varcurb
                                         if(printable.opt)print("backward Model Fit Error!?")
                                         next
                                       }
                                     }

                                     varcand<-vect[[mod_id]]$varcur
                                     if(cluster)
                                     {
                                       waiccand<-res.par[[mod_id]]$waic
                                       mlikcand<-res.par[[mod_id]]$mlik
                                     }else if(exists("statistics1"))
                                     {
                                       iidd<-bittodec(varcand)+1
                                       waiccand<-statistics1[iidd,2]
                                       mlikcand<-statistics1[iidd,1]
                                     }else if(exists("hashStat"))
                                     {
                                       iidd<-paste(varcand,collapse = "")
                                       waiccand<-values(hashStat[iidd])[2]
                                       mlikcand<-values(hashStat[iidd])[1]
                                     }

                                     if(objective==0)
                                     {
                                       objcand<-waiccand
                                       objcur<-waiccur
                                       objglob<-waicglob
                                     }else
                                     {
                                       objcand<- -mlikcand
                                       objcur<-  -mlikcur
                                       objglob<- -mlikglob
                                     }


                                     if(objcand<=objcur|| mod_id ==1)
                                     {
                                       if(printable.opt)print(paste("backward accept with ", objcand))
                                       objcur<-objcand
                                       varcurb<-varcand
                                       waiccur<-waiccand
                                       varcur<-varcand
                                       mlikcur<-mlikcand

                                       if(objcur<objglob)
                                       {
                                         objglob<-objcur
                                         waicglob<-waiccand
                                         varglob<-varcand
                                         mlikglob<-mlikcand
                                         i