rep <- 200
library(OpenMx)
library(doParallel)

cat( toLatex(sessionInfo(), locale = FALSE) , file="rsessioninfo.txt", sep = "\n" )

sim <- function(beta1=.14, beta2=.14, samp=100){
b1 <- beta1;b2 <- beta2; n <- samp
s2x <- 1; s2m <- 1; s2y <- 1
manifests = c("x", "m", "y")

popModel <- mxModel(
  name = "populationModel",
  type = "RAM",
  manifestVars = manifests,
  mxPath(
    from = 'one',
    to = manifests,
    free = TRUE,
    values = 0,
    labels = c("mu_x", "int_m", "int_y")
  ),
  mxPath(
    from = "x",
    to = "m",
    arrows = 1,
    free = TRUE,
    values = b1,
    labels = "b1"
  ),
  mxPath(
    from = "m",
    to = "y",
    arrows = 1,
    free = TRUE,
    values = b2,
    labels = "b2"
  ),
  mxPath(
    from = manifests,
    arrows = 2,
    free = TRUE,
    values = c(s2x, s2m, s2y),
    labels = c("s2x", "s2m", "s2y")
  )
)
popModel <- mxOption(popModel, "Standard Errors"  , "No")
## Generating Data and Retrun the model that genberates Data
sim_df <-
  mxGenerateData(
    model = popModel,
    nrows = n,
    returnModel = FALSE,
    use.miss = FALSE
  )

medModel <-
  mxModel(model = popModel,
    name = "Unconstrained Model",
    mxData(sim_df, type = "raw"))

## Constrained Model
medModelconst <- mxModel(
  model = medModel,
  name = "constrained Model",
  mxAlgebra(b1 * b2, name = "ind"),
  mxConstraint(ind == 0, name = "ind0")
)

medModelconstFit <-
  mxTryHard(medModelconst,
    extraTries = 10,
    checkHess = FALSE)

medModelFit <- mxRun(medModel)


cbind(
  pop = omxGetParameters(popModel),
  uncosnt = omxGetParameters(medModelFit),
  const = omxGetParameters(medModelconstFit)
)

}

set.seed(12387, "L'Ecuyer")
mc <- detectCores() - 1
cl <- makeCluster(mc)
registerDoParallel(cl)

it <- 1:rep
resSim <-
  foreach(
    x = it,
    .combine = "list",
    .packages = c("OpenMx"),
    .inorder = FALSE,
    .errorhandling = "stop",
    .verbose = TRUE
  ) %dopar% sim()

stopCluster(cl)

write.csv(resSim,"results.csv")