Cacher is the code snippet organizer for pro developers

library(bnlearn)
getR <- function(k, S){
  #Simulate a valid covariance matrix from an inverse Wishart
  tryCatch(#This uses the inverse of a wishart sim, non-zero prob of compuational singlularity,
    riwish(k, S),
    error = function(e) getR(k, S)
  )
}

simData <- function(k, n, var.names){
  #Simulate multivariate binomial with covariance
  require(MCMCpack)
  
  S <- toeplitz((k:1))
  resample <- sample(k, k)
  S <- S[resample, resample]
  R <- getR(k, S)
  sim.data <- tryCatch(
    as.data.frame(matrix(rnorm(n * k), ncol = k) %*% chol(R) > 0),
    error = function(e) simData(k, n, var.names)
  )
  for(i in 1:ncol(sim.data)){
    sim.data[, i] <- as.factor(sim.data[, i])
    levels(sim.data[, i]) <- c("LOW", "HIGH")
  }  
  names(sim.data) <- var.names
  sim.data
}
k <- 5 # number of  nodes
#Simulate a network structure
node.names <- paste("S", 1:k, sep = "")
node.net <-  random.graph(nodes = node.names, method = "melancon", num = 1, burn.in = 10^5, every = 100)[[1]]

#Fit parameters to the network
base.data <- simData(k, 10, nodes(node.net))
sim.data <- rbn(node.net, 100, base.data)
fitted.net <- bn.fit(node.net, sim.data, method="bayes")

#Recursively simulate states for all nodes in the network
simState <- function(x, parent.states){
  #Simulate a state for a given node given the parents
  parent.states
  sample.probs <- eval(parse(text = paste(c("probs.list[[x]][,", paste(paste("'", parent.states, "'", sep = ""), collapse = ","), "]"), collapse = "")))
  sample(c("LOW", "HIGH"), 1, prob = sample.probs)
}
simLineage <- function(time.t.states, node){
  #Recursively simulate the states of each node.  Given a node, simulate a state given its parents.
  #If the parents do not yet have a state, recursively apply the function to each parent
  node.parents <- parents(node.net, node)
  not.yet.simmed <- node.parents[sapply(node.parents, function(p) is.na(time.t.states[p]))]
  if(length(not.yet.simmed) > 0){
    for(p in not.yet.simmed){
      time.t.states <- simLineage(time.t.states, p)
    }
  }
  
  parent.states <- time.t.states[node.parents]
  time.t.states[node] <- simState(node, parent.states)
  time.t.states
}


#Make list of roots (nodes with no parents) and youngest (nodes with no children)
parents.list <- sapply(nodes(node.net), function(node) parents(node.net, node))
roots <- nodes(node.net)[sapply(parents.list, function(item) length(item) == 0)]
children.list <- sapply(nodes(node.net), function(node) children(node.net, node))
youngest <- nodes(node.net)[sapply(children.list, function(item) length(item) == 0) ]
probs.list <- lapply(nodes(node.net), function(node) fitted.net[[node]]$prob)
names(probs.list) <- nodes(node.net)
time.t.states <- rep(NA, nnodes(node.net))
names(time.t.states) <- nodes(node.net)

#Roots will be the base case in the recursion.  Simulate their states first.
time.t.states[roots] <- sapply(probs.list[roots], function(prob.vals) sample(c("LOW", "HIGH"), 1, prob=prob.vals))

#Starting with the youngest, apply the recursion such that it propagates back through the network until the roots are reached.
for(y in youngest){
  time.t.states <- simLineage(time.t.states, y)
}
time.t.states