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

## ----sweep--------------------------------------------------------------------
sim <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 11)
Y <- sim$dat
groups <- rep(1:3, each = 6)

Q0 <- matrix(-1L, ncol(Y), 2)
for (k in 1:2) Q0[which(groups == k)[1:2], k] <- 1L
Q0[1:8, ]

p <- pefa(Q0, Y, Kmin = 2, Kmax = 7, verbose = FALSE)
p

## ----object-------------------------------------------------------------------
names(p)
class(p)
names(p$settings)

## ----matrices-----------------------------------------------------------------
names(p$loadings)
names(p$pips)
round(p$loadings[["3"]], 2)
round(p$pips[["3"]][1:6, ], 2)

## ----sweep-table--------------------------------------------------------------
p$sweep

## ----transitions--------------------------------------------------------------
p$transitions

## ----persistence--------------------------------------------------------------
names(p$persistence)
round(p$persistence$phi, 3)
round(p$persistence$rmsd, 3)
p$persistence$collision

## ----k2-----------------------------------------------------------------------
round(p$loadings[["2"]][c(1, 7, 8, 13, 17, 18), ], 2)

## ----ssl----------------------------------------------------------------------
ssl(p)

## ----k7-----------------------------------------------------------------------
round(p$loadings[["7"]][13:18, 3:5], 2)

## ----reading------------------------------------------------------------------
profile_specs <- data.frame(
  profile = c(".85/r1", ".80/r2", ".70/r3"),
  phi_cut = c(.85, .80, .70), r = 1:3,
  practical_default = c(FALSE, TRUE, FALSE)
)
profile_sources <- 2:4

read_sweep <- function(x, phi_cut, r, sources) {
  K <- as.integer(x$sweep$K)
  converged <- x$sweep$converged

  read_count <- function(score) {
    usable <- length(K) >= 2L && all(is.finite(K)) && all(diff(K) == 1L) &&
      length(score) == length(K) && all(is.finite(score)) &&
      is.logical(converged) && all(converged %in% TRUE)
    if (!usable) return(list(usable = FALSE, Khat = NA_integer_))
    gain <- diff(score)
    gain_max <- max(gain)
    if (!(gain_max > 0)) return(list(usable = TRUE, Khat = NA_integer_))
    peak <- max(which(gain == gain_max))
    hit <- which(seq_along(gain) > peak & gain < .20 * gain_max)
    list(usable = TRUE,
         Khat = if (length(hit)) K[hit[1L]] else NA_integer_)
  }

  ELBO <- read_count(x$sweep$ELBO)
  BIC <- read_count(-x$sweep$BIC)
  C20 <- sort(unique(c(ELBO$Khat, BIC$Khat)[
    is.finite(c(ELBO$Khat, BIC$Khat))]))

  phi <- x$persistence$phi
  collision <- x$persistence$collision
  conv <- setNames(converged, K)
  state <- vapply(sources, function(k) {
    targets <- k + seq_len(r)
    edge <- vapply(targets, function(target) {
      from <- as.character(k); to <- as.character(target)
      if (!from %in% rownames(phi) || !to %in% colnames(phi) ||
          !isTRUE(conv[[from]]) || !isTRUE(conv[[to]]) ||
          is.na(collision[from, to])) return(NA_integer_)
      if (isTRUE(collision[from, to])) return(0L)
      value <- phi[from, to]
      if (!is.finite(value)) NA_integer_ else as.integer(value >= phi_cut)
    }, integer(1))
    if (any(edge %in% 0L)) 0L else
      if (length(edge) && all(edge %in% 1L)) 1L else NA_integer_
  }, integer(1))

  persistent <- which(state == 1L)
  if (!length(state)) {
    resolution <- "unresolved"
    Kp <- NA_integer_
  } else if (!length(persistent)) {
    resolution <- if (all(state %in% 0L)) "nonpersistent" else "unresolved"
    Kp <- NA_integer_
  } else {
    top <- max(persistent)
    higher <- if (top < length(state)) state[(top + 1L):length(state)] else 0L
    resolution <- if (all(higher %in% 0L)) "persistent" else "unresolved"
    Kp <- if (resolution == "persistent") sources[top] else NA_integer_
  }

  hit_ELBO <- is.finite(Kp) && is.finite(ELBO$Khat) && ELBO$Khat == Kp
  hit_BIC <- is.finite(Kp) && is.finite(BIC$Khat) && BIC$Khat == Kp
  support <- if (hit_ELBO && hit_BIC) "both" else if (hit_ELBO) "ELBO only" else
    if (hit_BIC) "BIC only" else "none"
  any_count_usable <- ELBO$usable || BIC$usable
  layer <- if (resolution == "persistent" && Kp %in% C20) "L1" else
    if (resolution == "persistent" && any_count_usable) "L2" else
    if (resolution == "nonpersistent") "L3" else
      "unclassified"
  gap <- if (resolution == "persistent" && length(C20))
    min(abs(C20 - Kp)) else NA_real_

  list(Khat_ELBO = ELBO$Khat, Khat_BIC = BIC$Khat, C20 = C20,
       count_usable = c(ELBO = ELBO$usable, BIC = BIC$usable),
       support = support, Kp = Kp, source_state = setNames(state, sources),
       resolution = resolution, layer = layer, g_CP = gap)
}

profile_row <- function(x, spec) {
  z <- read_sweep(x, spec$phi_cut, spec$r, profile_sources)
  data.frame(profile = spec$profile, practical_default = spec$practical_default,
             Khat_ELBO = z$Khat_ELBO, Khat_BIC = z$Khat_BIC,
             C20 = if (length(z$C20)) paste0("{", paste(z$C20,
               collapse = ","), "}") else "{}",
             support = z$support, Kp = z$Kp, layer = z$layer,
             g_CP = z$g_CP, row.names = NULL)
}

profile_results <- do.call(rbind, lapply(seq_len(nrow(profile_specs)),
  function(i) profile_row(p, profile_specs[i, ])))
profile_results

## A declared source at the upper edge cannot persist without its endpoint.
short_window <- read_sweep(p, phi_cut = .85, r = 1L, sources = 7L)
data.frame(profile = ".85/r1", declared_source = 7L,
           required_endpoint = 8L, resolution = short_window$resolution,
           layer = short_window$layer)

## ----summary------------------------------------------------------------------
summary(p)

## ----plots, fig.width=7, fig.height=4.5---------------------------------------
plot(p, type = "objective", criterion = "ELBO")
plot(p, type = "gain", criterion = "BIC", pct = TRUE)
plot(p, type = "fit")