---
title: "Multistart Penalized Estimation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multistart Penalized Estimation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```

```{r setup}
library(plavaan)
library(lavaan)
data(PoliticalDemocracy)
```

## Why multistart?

Penalized objectives such as `l0a` and `alf` are non-convex near zero. When combined with a generic optimizer like `nlminb()`, the risk of settling in a local optimum is real. Running the optimization from multiple starting values helps avoid this problem: if different starts converge to the same solution, you can be more confident it is (locally) optimal; if they diverge, the multistart wrapper selects the best among them.

## Single-start baseline

Start with the longitudinal invariance example from [penalized_est()]. We fit a two-wave CFA and penalize differences in loadings and intercepts across time:

```{r 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 with random jittering

`penalized_est_multistart()` wraps `penalized_est()` and tries several starting vectors. The first start is always lavaan's default (unperturbed), so multistart can never do worse than the single-start baseline:

```{r 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")
```

The `multistart` attribute is a data frame with one row per start, showing the final penalized objective (`objective`) and whether it converged. Rows are sorted by ascending objective.

## Comparing local solutions

When starts lead to different parameter estimates, it can be useful to inspect all solutions:

```{r 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
```

## Using custom starting values

Instead of random jittering, you can supply your own starting vectors via `starts`:

```{r 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")
```

## Reproducibility

The function does not call `set.seed()` internally. Set your own seed before calling for reproducibility:

```{r 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 execution

Multistart runs sequentially for simplicity. For parallel execution, call `penalized_est(start = ...)` directly with your preferred backend (`future.apply`, `parallel`, etc.):

```{r 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, ]
    )
})
```