## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(plavaan)
library(lavaan)
data(PoliticalDemocracy)

## ----single-start-------------------------------------------------------------
model <- "
  dem60 =~ y1 + y2 + y3 + y4
  dem65 =~ y5 + y6 + y7 + y8
  dem60 ~~ dem65
  dem60 ~~ 1 * dem60
  dem65 ~~ NA * dem65
  dem60 ~ 0
  dem65 ~ NA * 1
  y1 ~~ y5
  y2 ~~ y6
  y3 ~~ y7
  y4 ~~ y8
"

fit_un <- cfa(model, data = PoliticalDemocracy, std.lv = TRUE,
              meanstructure = TRUE, do.fit = FALSE)

pt <- parTable(fit_un)
load_60 <- pt$free[pt$op == "=~" & pt$lhs == "dem60"]
load_65 <- pt$free[pt$op == "=~" & pt$lhs == "dem65"]
int_60  <- pt$free[pt$op == "~1" & pt$lhs %in% c("y1", "y2", "y3", "y4")]
int_65  <- pt$free[pt$op == "~1" & pt$lhs %in% c("y5", "y6", "y7", "y8")]

fit_single <- penalized_est(
    x = fit_un,
    w = 0.03,
    pen_diff_id = list(
        loadings   = rbind(load_60, load_65),
        intercepts = rbind(int_60, int_65)
    ),
    pen_fn = "l0a"
)

fit_single@optim$fx

## ----multistart---------------------------------------------------------------
set.seed(3821)

fit_multi <- penalized_est_multistart(
    x = fit_un,
    w = 0.03,
    pen_diff_id = list(
        loadings   = rbind(load_60, load_65),
        intercepts = rbind(int_60, int_65)
    ),
    pen_fn = "l0a",
    n_starts = 10,
    verbose = TRUE
)

# Inspect the spread of objective values across starts
attr(fit_multi, "multistart")

## ----keep-all-----------------------------------------------------------------
set.seed(3821)

fit_all <- penalized_est_multistart(
    x = fit_un,
    w = 0.03,
    pen_diff_id = list(
        loadings   = rbind(load_60, load_65),
        intercepts = rbind(int_60, int_65)
    ),
    pen_fn = "l0a",
    n_starts = 5,
    keep_all = TRUE
)

# Extract loadings from the top-3 solutions for comparison
all_fits <- attr(fit_all, "all_fits")
ms_table <- attr(fit_all, "multistart")[order(attr(fit_all, "multistart")$objective), ]

loadings_top3 <- lapply(head(all_fits[ms_table$start_id[1:3]], 3), function(f) {
    coefs <- coef(f)[grep("^dem60=~", names(coef(f)))]
    round(coefs, 4)
})
names(loadings_top3) <- paste0("Start ", ms_table$start_id[1:3])
loadings_top3

## ----custom-starts------------------------------------------------------------
base <- lavaan::lav_export_estimation(fit_un)$starting_values
my_starts <- rbind(
    base,                              # unperturbed
    base + runif(length(base), -0.1, 0.1)  # perturbed
)

fit_custom <- penalized_est_multistart(
    x = fit_un,
    w = 0.03,
    pen_diff_id = list(
        loadings   = rbind(load_60, load_65),
        intercepts = rbind(int_60, int_65)
    ),
    starts = my_starts
)

attr(fit_custom, "multistart")

## ----seed---------------------------------------------------------------------
set.seed(7293)
fit_repro <- penalized_est_multistart(
    x = fit_un, w = 0.03,
    pen_diff_id = list(loadings = rbind(load_60, load_65)),
    n_starts = 5
)

## ----parallel-pattern, eval = FALSE-------------------------------------------
# library(future.apply)
# plan(multicore)
# 
# # Note: random_start() is an internal helper without API stability guarantees.
# starts <- plavaan:::random_start(fit_un, n = 10)[-1, , drop = FALSE]  # skip base
# fits <- future_lapply(seq_len(nrow(starts)), function(i) {
#     penalized_est(
#         x = fit_un, w = 0.03,
#         pen_diff_id = list(loadings = rbind(load_60, load_65)),
#         start = starts[i, ]
#     )
# })