n <- c(50, 250, 1000)

for(i in 1:length(n)) dir.create(paste0("n", n[i]))

nRep <- 1000
currentDir <- getwd()

set.seed(123321)
for(i in 1:length(n)) {
  for(j in 1:nRep) {
    x <- rnorm(250, 0, 1)
	m1 <- 0.2 * x + rnorm(250, 0, sqrt(0.96))
	m2 <- 0.2 * x + rnorm(250, 0, sqrt(0.96))
	y <- 0.2 * x + 0.2 * m1 + 0.2 * m2 + rnorm(250, 0, sqrt(0.84)) 
	dat <- data.frame(x, m1, m2, y)
    filename <- paste0("dat", j, ".csv")
    targetFile <- paste0(currentDir, "/n", n[i], "/", filename)
    write.table(dat, file = targetFile, row.names = FALSE, col.names = FALSE, sep = ",")
  }
}

template <- "
  TITLE: Example on Mutliple-Mediation Model;
  DATA: 
    FILE IS subFile;
  VARIABLE: 
    NAMES ARE x m1 m2 y;
  MODEL:  
    y ON x (c) 
         m1 (b1)
		 m2 (b2);
    m1 ON x (a1);
	m2 ON x (a2);
  MODEL CONSTRAINT:
    new(ab1, ab2, indirect);
    ab1 = a1 * b1;
	ab2 = a2 * b2;
	indirect = ab1 + ab2;
"
for(i in 1:length(n)) {
  for(j in 1:nRep) {
    datname <- paste0("dat", j, ".csv")
    tempscript <- gsub("subFile", datname, template)
    filename <- paste0("script", j, ".inp")
    targetFile <- paste0(currentDir, "/n", n[i], "/", filename)
    write(tempscript, targetFile)
  }
}

for(i in 1:length(n)) {
  targetDir <- paste0(currentDir, "/n", n[i], "/")
  setwd(targetDir)
  for(j in 1:nRep) {
    filename <- paste0("script", j, ".inp")
    outname <- paste0("output", j, ".out")
    shell(paste0("mplus ", filename, " ", outname))
  }
}
setwd(currentDir)

result <- NULL
for(i in 1:length(n)) {
  for(j in 1:nRep) {
    targetFile <- paste0(currentDir, "/n", n[i], "/", "output", j, ".out")
    tempoutput <- readLines(targetFile)
    if(length(grep("THE MODEL ESTIMATION TERMINATED NORMALLY", tempoutput)) > 0) {
      targetline <- grep("INDIRECT", tempoutput)
      textline <- tempoutput[targetline]
      textsplit <- strsplit(textline, " ")
      textresult <- setdiff(textsplit[[1]], "")
      result <- rbind(result, c(n[i], as.numeric(textresult[2:5])))
    }
  }
}
colnames(result) <- c("n", "est", "se", "z", "p")
aggregate(est ~ n, data = result, FUN = mean)