manifest_names <- names(data_wide[ , 2:ncol(data_wide)])
return(manifest_names)
}

#create nonlinear latent growth curve model
create_logistic_growth_model <- function(data_wide, model_name, starting_values) {

manifest_vars <- generate_manifest_var_names(data = data_wide)
latent_vars <- c('diff_fixed', 'beta_fixed', 'gamma_fixed')
measurement_days <- as.numeric(str_extract(string = names(data_wide[ 2:ncol(data_wide)]), pattern = '[^_]*$'))

model <- mxModel(model = model_name,
type = 'RAM', independent = T,
mxData(observed = data_wide, type = 'raw'),

manifestVars = manifest_vars,
latentVars = latent_vars,

#Manifest means
manMeans = mxPath(from = 'one', to = manifest_vars, free = FALSE, arrows = 1),

#Latent variable means (linear parameters). Note that the nonlinear parameters of beta and gamma do not have estimated means
mxPath(from = 'one', to = c('diff_fixed'),free = TRUE, arrows = 1, values = starting_values$diff_fixed,
labels = 'diff_fixed', lbound = 0, ubound = 2),

#Residual variances; by using one label, they are assumed to all be equal (homogeneity of variance)
mxPath(from = manifest_vars,
arrows=2, free=TRUE, values = starting_values$epsilon, labels='epsilon', lbound = 1e-3),

#Latent variable covariances and variances
mxPath(from = latent_vars,
connect='unique.pairs', arrows=2,
#aa(diff_rand), ab(cov_diff_beta), ac(cov_diff_gamma), bb(beta_rand), bc(var_beta_gamma), cc(gamma_rand)
free = c(TRUE,
FALSE, FALSE,
TRUE, FALSE, TRUE),
values=c(starting_values$diff_rand,
NA,NA,
starting_values$beta_rand, NA, starting_values$gamma_rand),
labels=c('diff_rand',
'NA(cov_diff_beta)','NA(cov_diff_gamma)',
'beta_rand', 'NA(var_beta_gamma)', 'gamma_rand'),
lbound = c(1e-3,
NA, NA,
1, NA, 1)),

#Functional constraints
mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE,
values = starting_values$diff_fixed, labels = 'diff_fixed', name = 'd', lbound = 0, ubound = 2),
mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE,
values = starting_values$beta_fixed, labels = 'beta_fixed', name = 'b', lbound = 0, ubound = 360),
mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE,
values = starting_values$gamma_fixed, labels = 'gamma_fixed', name = 'g', lbound = 0, ubound = 360),
mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = FALSE, values = measurement_days, name = 'time'),

mxAlgebra(expression = 1/(1 + exp((beta_fixed - time)/gamma_fixed)), name="Dl"),
mxAlgebra(expression = -(diff_fixed*(exp((beta_fixed - time)/gamma_fixed) *
(1/gamma_fixed))/(1 + exp((beta_fixed - time)/gamma_fixed))^2), name = 'Bl'),
mxAlgebra(expression = diff_fixed * (exp((beta_fixed - time)/gamma_fixed) * ((beta_fixed - time)/gamma_fixed^2))/
(1 + exp((beta_fixed - time)/gamma_fixed))^2, name = 'Gl'),

#Factor loadings; all fixed and, importantly, constrained to equal their partial derivatives
mxPath(from = 'diff', to = manifest_vars, arrows=1, free=FALSE,
labels = c("Dl[1,1]", "Dl[2,1]", "Dl[3,1]", "Dl[4,1]", "Dl[5,1]")),
mxPath(from='beta', to = manifest_vars, arrows=1, free=FALSE,
labels = c("Bl[1,1]", "Bl[2,1]", "Bl[3,1]", "Bl[4,1]", "Bl[5,1]")),
mxPath(from='gamma', to = manifest_vars, arrows=1, free=FALSE,
labels = c("Gl[1,1]", "Gl[2,1]", "Gl[3,1]", "Gl[4,1]", "Gl[5,1]")),

mxFitFunctionML()
)

return(model)
}

#changes starting values of model created from create_logistic_growth_model
update_starting_values <- function(model_number, model, starting_values){

param_name_order <- names(omxGetParameters(model))

#set model with new starting values
model_name <- paste("Iteration", model_number, sep="")
new_model <- omxSetParameters(model = model,
name = model_name,
labels = param_name_order,
values = as.numeric(setcolorder(x = starting_values[model_number, ],
neworder = param_name_order)))

#return model with new starting values
return(new_model)
}

#iterates update_starting_values to create a parent model where each submodel has different starting values
assemble_model_starts <- function(model, starting_values) {

parent_model <- mxModel(name = "parent_model")
sub_models <- lapply(X = 1:nrow(starting_values), FUN = update_starting_values, model = model, starting_values = starting_values)
merged_model <- mxModel(parent_model, sub_models)

return(merged_model)
}

#extract necessary information about fit from each submodel
extract_fit_info <- function(model){

#read in starting values and raw data sets
starting_values <- read_csv('starting_values.csv')
data_wide <- read_csv('data_wide.csv')

#create nonlinear latent growth curve model
latent_growth_model <- create_logistic_growth_model(data_wide = data_wide,
model_name = 'lg_nonlinear', starting_values = starting_values[1, ])

#run submodel creation and model fitting in parallel mode
sfInit(parallel = TRUE, cpus = detectCores() - 1) #save 1 core to prevent computer from crashing
sfLibrary(OpenMx)

#create submodels that each have a different set of starting values
merged_models <- assemble_model_starts(model = latent_growth_model, starting_values = starting_values)

#run all models
output <- mxRun(merged_models, unsafe = T)

resultsFit <- data.table(t(omxSapply(output$submodels, fitStats)))
names(resultsFit) <- c("fit", "status_code", "num_iterations")
sum(resultsFit$status_code == 0) #0

latent_growth_model_ideal <- create_logistic_growth_model(data_wide = data_wide,
model_name = 'lg_nonlinear', starting_values = ideal_starting_values)

first_run <- mxRun(latent_growth_model_ideal) #status code 6

mxRun(first_run) #status code 5