R/nmf.ffb.R

In nmfkc: Non-Negative Matrix Factorization with Kernel Covariates

# nmf.ffb.R — NMF-FFB (formerly NMF-SEM) canonical engines + generic DOT
# Canonical:  nmf.ffb, nmf.ffb.inference, nmf.ffb.cv, nmf.ffb.split, nmf.ffb.DOT
#             (the deprecated nmf.sem* aliases were removed in 0.9.8).
# Also hosts: nmfkc.DOT / plot.nmfkc.DOT (shared DOT utilities).

#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#' @title NMF-FFB: non-negative matrix factorization with latent feedback
#'
#' @description
#' Fits the NMF-FFB model
#' \deqn{
#'   Y_1 = X B + E, \qquad B = \Theta_1 Y_1 + \Theta_2 Y_2 + U,
#' }
#' in which a non-negative basis \eqn{X} (\eqn{P_1 \times Q}) generates the
#' endogenous block \eqn{Y_1} from latent scores \eqn{B} that are driven by
#' the exogenous block \eqn{Y_2} (feed-forward, \eqn{\Theta_2}) and fed back
#' by \eqn{Y_1} itself (latent feedback, \eqn{\Theta_1}).  The function
#' returns the estimated basis, the structural coefficient matrices and the
#' implied equilibrium (input-output) mapping.
#'
#' At equilibrium, the model can be written as
#' \deqn{
#'   Y_1 \approx (I - X \Theta_1)^{-1} X \Theta_2 Y_2
#'   \equiv M_{\mathrm{model}} Y_2,
#' }
#' where \eqn{M_{\mathrm{model}} = (I - X \Theta_1)^{-1} X \Theta_2} is a
#' Leontief-type cumulative-effect operator in latent space.
#'
#' Internally, the latent feedback and exogenous loading matrices are stored as
#' \code{C1} and \code{C2}, corresponding to \eqn{\Theta_1} and \eqn{\Theta_2},
#' respectively.
#'
#' Two estimators are available through \code{method}:
#' \describe{
#'   \item{\code{"fiml"} (default)}{A two-stage likelihood-based estimator.
#'     \strong{Stage 1} estimates the basis with the feed-forward fit
#'     \code{\link{nmfkc}(Y1, A = Y2)} (\code{X.init}, \code{X.L2.ortho},
#'     \code{epsilon}, \code{maxit}, \code{seed} are forwarded), unless a
#'     basis is supplied through \code{X}.  Columns are normalized to unit
#'     sum.  \strong{Stage 2} holds \eqn{\hat X} fixed and fits the Gaussian
#'     working model \eqn{U \sim N(0, \Phi)}, \eqn{E \sim N(0,
#'     \mathrm{diag}(\psi))} by full-information maximum likelihood
#'     (L-BFGS-B with an analytic gradient) under the non-negativity of
#'     \eqn{\Theta_1, \Theta_2} and an exclusion restriction on
#'     \eqn{\Theta_1}.  It fits (a) the feed-forward null \eqn{\Theta_1 = 0}
#'     (a non-negative MIMIC factor model with correlated factors), (b) the
#'     unpenalized feedback model on the admitted entries, and (c) an L1
#'     path over \code{C1.L1.path}.  The penalized problem is non-convex, so
#'     every point of the path is fitted from several starting points
#'     (\code{starts}); every distinct support proposed by any (penalty,
#'     start) pair is re-estimated without penalty (itself from three
#'     starts, keeping the best log-likelihood) and BIC is minimized over
#'     all distinct candidates together with the null and the unpenalized
#'     model.  The support with the smallest BIC is the \emph{selected}
#'     model, reported in \code{C1}, \code{C2}, \code{Phi}, \code{psi}.
#'     \code{C1.L1} and \code{C2.L1} are not used by this estimator.  The likelihood-ratio statistics against the null are
#'     returned without p-values: \eqn{\Theta_1 \ge 0} puts the null on the
#'     boundary of the parameter space and the BIC refit is a post-selection
#'     statistic, so a chi-square reference is not valid.  Their calibration
#'     by parametric bootstrap is done by \code{\link{nmf.ffb.inference}}.}
#'   \item{\code{"mu"}}{The legacy estimator: joint multiplicative updates of
#'     \eqn{X, \Theta_1, \Theta_2} for the structural-form squared error
#'     \eqn{\lVert Y_1 - X(\Theta_1 Y_1 + \Theta_2 Y_2)\rVert_F^2} with an
#'     orthogonality penalty on \eqn{X} and L1 penalties \code{C1.L1},
#'     \code{C2.L1}.  Kept, bit-identical, so that analyses published with it
#'     reproduce.  Because the structural and reduced forms have the same
#'     fit once \eqn{X} is free, this estimator cannot separate
#'     \eqn{\Theta_1} from \eqn{\Theta_2}; prefer \code{"fiml"}.}
#' }
#'
#' @section Exclusion restriction:
#' Feedback is identified only through exclusion restrictions: an outcome
#' may not feed back into a factor on which it has more than a negligible
#' loading, because such an entry is nearly equivalent to a change of the
#' outcome's loading and uniqueness.  With
#' \code{C1.restriction = "union"} (default) entry \eqn{(q, i)} of \eqn{\Theta_1} is
#' excluded if \eqn{q = \arg\max_{q'} X_{i q'}} (the dominant factor) or
#' \eqn{X_{iq} \ge} \code{C1.restriction.threshold}.  Both halves are needed:
#' blocking the dominant factor alone would leave an outcome with a
#' substantial second loading free to feed that factor, and blocking only the
#' factors above the threshold would leave an outcome whose largest loading is
#' below the threshold free to feed its own.
#' \code{"none"} frees every entry (not recommended: the model is then
#' identified only through the non-negativity and the covariance structure).
#' A user-supplied \eqn{Q \times P_1} 0/1 matrix is used as given.
#'
#' The restriction is derived from the estimated basis and therefore from the same
#' \eqn{Y_1} that is subsequently tested; see the \emph{Calibration} section
#' of \code{\link{nmf.ffb.inference}} for what this implies.
#'
#' @section Lifecycle:
#' \code{method = "fiml"} became the default in version 0.9.8, as did
#' \code{C1.restriction = "union"} (earlier fiml fits blocked the dominant
#' factor only), and
#' \code{\link{nmf.ffb.inference}} gained the \code{calibration} argument.
#' \code{method = "mu"} is the legacy estimator, kept for the reproducibility
#' of published analyses; it will be deprecated in a later release.
#'
#' @param Y1 A non-negative numeric matrix of endogenous variables with
#'   \strong{rows = variables (P1), columns = samples (N)}.
#' @param Y2 A non-negative numeric matrix of exogenous variables with
#'   \strong{rows = variables (P2), columns = samples (N)}.
#'   Must satisfy \code{ncol(Y1) == ncol(Y2)}.
#' @param rank Integer; number of latent factors \eqn{Q}. If \code{NULL},
#'   \eqn{Q} is taken from a hidden argument in \code{...} or defaults to
#'   \code{nrow(Y2)}.
#' @param X.init Initialization strategy for the basis matrix
#'   \code{X} (\eqn{P_1 \times Q}).  One of:
#'   \itemize{
#'     \item \code{"nndsvd"} (default): Non-negative Double SVD with
#'       additive randomness (NNDSVDar; Boutsidis & Gallopoulos 2008),
#'       computed internally via \code{.nndsvdar(Y1, Q)}.  Requires
#'       \eqn{Q \le \min(P_1, N)} (over-rank case falls back to
#'       \code{"runif"}).  Uses a full SVD of \eqn{Y_1}, so for very
#'       large \eqn{Y_1} consider switching to \code{"kmeans"} to
#'       avoid SVD memory / compute cost.
#'     \item \code{"kmeans"}: k-means on the columns of \eqn{Y_1}
#'       (samples clustered into \eqn{Q} groups); the transposed
#'       cluster centers become \eqn{X}.  Scales well for large
#'       \eqn{Y_1}; this is the default of \code{\link{nmfkc}}.
#'     \item \code{"kmeansar"}: \code{"kmeans"} followed by filling
#'       zero entries of \eqn{X} with \eqn{\mathrm{Uniform}(0,
#'       \bar Y_1 / 100)} (NNDSVDar-style additive randomness to
#'       escape trivial stationary points).
#'     \item \code{"runif"}: Uniform random entries in \eqn{[0, 1]}.
#'     \item \code{"spa"}: the anchor rows the successive projection
#'       algorithm finds in \eqn{Y_1} (see \code{\link{nmfkc}}); a starting
#'       value only.  \code{method = "fiml"} only.
#'     \item A numeric \eqn{P_1 \times Q} matrix supplied by the user;
#'       negative entries are projected to 0.
#'     \item \code{NULL}: backward-compatible alias for \code{"nndsvd"}.
#'   }
#'   To impose anchor rows on the stage-1 basis instead, pass
#'   \code{X.anchor} (\code{"spa"} or one row of \eqn{Y_1} per basis) through
#'   \code{...}, as in \code{\link{nmfkc}}: stage 2 holds the basis fixed, so
#'   the zeros stay.  It needs \code{rank <= nrow(Y2)}, is not available with
#'   \code{method = "mu"} or a supplied \code{X}, and the rows used are
#'   returned as \code{X.anchor}.
#'   In all cases the result is column-normalized to \code{colSums(X) = 1}
#'   before iteration.  The menu mirrors \code{\link{nmfkc}}'s
#'   \code{X.init} option for consistency across the package.
#' @param X.L2.ortho L2 orthogonality penalty for \code{X}. This controls
#'   the penalty term \eqn{\lambda_X \lVert X^\top X - \mathrm{diag}(X^\top X)
#'   \rVert_F^2}. Default: \code{100}.
#' @param C1.L1 L1 sparsity penalty for \code{C1} (i.e., \eqn{\Theta_1}).
#'   Default: \code{1.0}.  \strong{Scale note}: this function adds
#'   \code{C1.L1} to the multiplicative denominator, whereas
#'   \code{\link{nmfkc}}, \code{\link{nmf.rrr}} and \code{\link{nmfkc.net}} add
#'   \code{C.L1 / 2}.  The same nominal value is therefore twice as strong
#'   here.  The difference is retained so that published \code{nmf.ffb} fits
#'   reproduce; halve the value to match the other models.
#' @param C2.L1 L1 sparsity penalty for \code{C2} (i.e., \eqn{\Theta_2}).
#'   Default: \code{0.1}.  Same scale note as \code{C1.L1}.
#' @param epsilon Relative convergence threshold for the objective function.
#'   Iterations stop when the relative change in \strong{reconstruction loss}
#'   falls below this value. Default: \code{1e-6}.
#'   \strong{Note}: the test is on the unpenalized loss (\code{objfunc}), not
#'   on the penalized objective the updates actually minimize
#'   (\code{objfunc.penalized}).  Every other optimizer in the package tests the
#'   penalized value, so with a large \code{X.L2.ortho} or \code{C*.L1} this
#'   function can stop while the quantity being optimized is still moving.
#'   Both traces are returned; compare them if the penalties are strong.
#' @param maxit Maximum number of iterations for the multiplicative updates.
#'   Default: \code{5000} (matches \code{\link{nmfkc}} and other MU
#'   functions in the package).
#' @param seed Random seed used to initialize \code{X}, \code{C1}, and \code{C2}.
#'   Default: \code{123}.  For \code{method = "fiml"} the seed only reaches
#'   the stage-1 \code{\link{nmfkc}} fit; the FIML stage is deterministic.
#' @param method \code{"fiml"} (default) or \code{"mu"}; see Description.
#' @param X Optional basis for \code{method = "fiml"}: a \eqn{P_1 \times Q}
#'   non-negative matrix, or an \code{\link{nmfkc}} / \code{nmf.ffb} object
#'   whose \code{$X} is used.  When supplied, stage 1 is skipped and
#'   \code{rank} is taken from \code{ncol(X)}.
#' @param C1.restriction Exclusion restriction on \eqn{\Theta_1} for
#'   \code{method = "fiml"}: \code{"union"} (default), \code{"none"}, or a
#'   \eqn{Q \times P_1} 0/1 matrix
#'   (1 = free).  See the section \emph{Exclusion restriction}.
#' @param C1.restriction.threshold Loading threshold for
#'   \code{C1.restriction = "union"}.  Default \code{0.05}.
#' @param Phi.restriction Covariance of the latent disturbance \eqn{U} for
#'   \code{method = "fiml"}: \code{"full"} (default; positive definite via
#'   Cholesky) or \code{"diag"}.
#' @param C1.L1.path Numeric vector of L1 penalties on \eqn{\Theta_1} defining
#'   the path (\code{method = "fiml"}).  Default \code{NULL}, meaning
#'   \code{N * c(0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)}.
#' @param select \code{"BIC"} (default): the support with the smallest BIC
#'   among all candidates proposed along the path is the selected model;
#'   \code{"none"}: the unpenalized feedback fit is returned as the selected
#'   model and the path is skipped.
#' @param ... Additional hidden arguments.  For \code{method = "fiml"}:
#'   \code{fiml.maxit} (L-BFGS-B iteration cap, default \code{3000}),
#'   \code{factr} (\code{optim} tolerance, default \code{1e3}), and
#'   \code{Q} (alias of \code{rank}).  For \code{method = "mu"} the
#'   following control the optional feedforward baseline (used both as an
#'   \eqn{X} warm-start and as the reference for \code{SC.map}, the
#'   input-output structural fidelity defined in Satoh (2025) §4.SC.map):
#'   \describe{
#'     \item{\code{nmfkc.baseline}}{Controls whether a feedforward
#'       \code{\link{nmfkc}}(Y1, A = Y2) fit is used as baseline.
#'       Possible values:
#'       \itemize{
#'         \item Default (not given) — \code{nmf.ffb} runs
#'           \code{\link{nmfkc}} \strong{internally} when \code{X.init}
#'           is a string method (\code{"nndsvd"}, \code{"kmeans"},
#'           \dots) or \code{NULL}, forwarding \code{X.init},
#'           \code{X.L2.ortho}, \code{epsilon}, \code{maxit},
#'           \code{seed}.  The fitted \eqn{X} of the baseline is then
#'           used as warm-start for the nmf.ffb MU iterations, and
#'           \code{SC.map} is computed.  This means
#'           \code{nmf.ffb(Y1, Y2, rank = Q)} runs end-to-end without
#'           a prior \code{nmfkc} call.
#'         \item \code{TRUE} — same as above, but force the internal
#'           \code{\link{nmfkc}} call even when \code{X.init} is a
#'           user-supplied matrix (the matrix is overridden).
#'         \item \code{FALSE} — opt out; no internal call,
#'           \code{SC.map = NA} (pre-v0.6.8 behavior).
#'         \item An \code{\link{nmfkc}} result (list with \code{$X}
#'           and \code{$C}) — use as the baseline directly (no
#'           internal call); also adopted as \code{X.init} when the
#'           latter is a string / NULL.
#'       }}
#'     \item{\code{M.simple}}{Optional \eqn{P_1 \times P_2} pre-computed
#'       baseline mapping.  Takes precedence over \code{nmfkc.baseline}
#'       for the SC.map calculation but does not affect warm-start.}
#'     \item{\code{Q}}{Backward-compat alias for \code{rank}.}
#'   }
#'
#' @return An object of class \code{c("nmf.ffb", "nmf")}, a list
#'   with components:
#'   \item{X}{Estimated basis matrix (\eqn{P_1 \times Q}).}
#'   \item{C1}{Estimated latent feedback matrix (\eqn{\Theta_1}, \eqn{Q \times P_1});
#'     for \code{"fiml"} the BIC-selected estimate.}
#'   \item{C2}{Estimated exogenous loading matrix (\eqn{\Theta_2}, \eqn{Q \times P_2}).}
#'   \item{XC1}{Feedback matrix \eqn{X \Theta_1}.}
#'   \item{XC2}{Direct-effect matrix \eqn{X \Theta_2}.}
#'   \item{XC1.radius}{Spectral radius \eqn{\rho(X \Theta_1)}.}
#'   \item{XC1.norm1}{Induced 1-norm \eqn{\lVert X \Theta_1 \rVert_{1,\mathrm{op}}}.}
#'   \item{Leontief.inv}{Leontief-type inverse \eqn{(I - X \Theta_1)^{-1}.}}
#'   \item{M.model}{Equilibrium mapping
#'     \eqn{M_{\mathrm{model}} = (I - X \Theta_1)^{-1} X \Theta_2}.}
#'   \item{amplification}{Latent amplification factor
#'     \eqn{\lVert M_{\mathrm{model}} \rVert_{1,\mathrm{op}} /
#'          \bigl\lVert X \Theta_2 \bigr\rVert_{1,\mathrm{op}}}; meaningful
#'     only when \code{XC1.radius > 0} (it equals 1 when no feedback is
#'     selected).}
#'   \item{amplification.bound}{Geometric-series upper bound
#'     \eqn{1 / (1 - \lVert X \Theta_1 \rVert_{1,\mathrm{op}})} if
#'     \eqn{\lVert X \Theta_1 \rVert_{1,\mathrm{op}} < 1}, otherwise \code{Inf}.}
#'   \item{rank}{Effective latent dimension used in the fit.}
#'   \item{SC.cov}{Correlation between sample and model-implied covariance
#'     (flattened) of \eqn{Y_1}.  See \emph{second-moment fidelity} in
#'     Satoh (2025).}
#'   \item{SC.map}{Correlation between the equilibrium operator
#'     \eqn{M_{\mathrm{model}}} and a feedforward baseline mapping
#'     \eqn{M_{\mathrm{simple}} = X_0 \Theta_0}, computed only when the
#'     baseline is supplied via \code{M.simple} or \code{nmfkc.baseline}
#'     in \code{...}; otherwise \code{NA}.  See \emph{input-output
#'     structural fidelity} in Satoh (2025).}
#'   \item{mae}{Mean absolute error between \eqn{Y_1} and its equilibrium
#'     prediction \eqn{\hat Y_1 = M_{\mathrm{model}} Y_2}.}
#'   \item{objfunc}{Vector of reconstruction losses per iteration
#'     (\code{"mu"}); \code{NULL} for \code{"fiml"}.}
#'   \item{objfunc.penalized}{Vector of penalized objective values per iteration
#'     (\code{"mu"}); \code{NULL} for \code{"fiml"}.}
#'   \item{iter, maxit, epsilon, converged}{Convergence bookkeeping.  For
#'     \code{"fiml"}: the number of objective evaluations of the selected
#'     fit, the L-BFGS-B cap (\code{fiml.maxit}), the stage-1 tolerance,
#'     and \code{optim()$convergence == 0}.}
#'   \item{method}{\code{"fiml"} or \code{"mu"}.}
#'
#'   The following are present for \code{method = "fiml"} only
#'   (\code{SC.cov} and \code{SC.map} are then \code{NULL}):
#'   \item{Phi, psi, loglik, npar}{Latent-disturbance covariance
#'     (\eqn{Q \times Q}), unique variances (length \eqn{P_1}),
#'     log-likelihood and parameter count of the selected model.}
#'   \item{null}{The feed-forward null: list \code{C2, Phi, psi, loglik,
#'     npar, M.model}.}
#'   \item{full}{The unpenalized feedback fit: list \code{C1, C2, Phi, psi,
#'     loglik, npar, XC1.radius}.}
#'   \item{path}{Data frame with one row per (\code{C1.L1.path}, start) of
#'     the L1 path (\code{C1.L1.path = 0} is the unpenalized fit re-estimated
#'     on its non-zero entries, \code{Inf} the null): \code{C1.L1.path, start,
#'     support_id, nnz, rho, loglik, BIC, MAE} (each after re-estimation on
#'     the proposed support),
#'     \code{pen.value} (the penalized objective reached by that start,
#'     smaller is better) and \code{duplicate} (\code{TRUE} when the same
#'     support was already proposed by an earlier row).}
#'   \item{candidates, supports, support.selected}{One row per distinct
#'     support (\code{support_id, nnz, rho, loglik, BIC, MAE,
#'     lambda1.first, start.first, selected}; the null and the model with
#'     every admitted entry free are always candidates, so
#'     \code{candidates} can hold a support that appears in no row of
#'     \code{path}); the supports themselves
#'     (list of logical \eqn{Q \times P_1} matrices indexed by
#'     \code{support_id}); and the id of the selected one.}
#'   \item{C1.free, C1.L1.path, C1.L1.selected, support}{The free-entry matrix used
#'     (\eqn{Q \times P_1} 0/1), the path, the smallest penalty at which
#'     the selected support was proposed (\code{0} for the unpenalized
#'     model, \code{Inf} for the null), the selected support (logical
#'     \eqn{Q \times P_1}) and the starts used.}
#'   \item{LR, LR.df}{Likelihood-ratio statistics \code{c(full = 2(l_full -
#'     l_null), selected = 2(l_sel - l_null))} and the naive degrees of
#'     freedom \code{c(full = sum(C1.restriction), selected = nnz)} (also stored as
#'     \code{attr(LR, "df")}).  No p-value is attached; see
#'     \code{\link{nmf.ffb.inference}}.}
#'   \item{BIC, AIC}{Named vectors \code{c(null, full, selected)}.}
#'   \item{call}{The matched call, from which
#'     \code{\link{nmf.ffb.inference}} inherits the design.}
#'
#' @examples
#' # Simple NMF-FFB with iris data (non-negative)
#' Y <- t(iris[, -5])
#' Y1 <- Y[1:2, ]  # Sepal
#' Y2 <- Y[3:4, ]  # Petal
#' result <- nmf.ffb(Y1, Y2, rank = 2)
#' result$LR          # feedback vs feed-forward null (no p-value here)
#' result$BIC
#' result$path
#'
#' # Legacy multiplicative-update estimator
#' result.mu <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500, method = "mu")
#' result.mu$mae
#'
#' @seealso \code{\link{nmf.ffb.inference}}, \code{\link{nmf.ffb.cv}},
#'   \code{\link{nmf.ffb.split}}, \code{\link{nmf.ffb.DOT}},
#'   \code{\link{summary.nmf.ffb}}
#' @references
#' Satoh, K. (2025). Applying non-negative matrix factorization with covariates
#'   to structural equation modeling for blind input-output analysis.
#'   arXiv:2512.18250. \url{https://arxiv.org/abs/2512.18250}
#' @export
nmf.ffb <- function(
    Y1, Y2,
    rank = NULL,
    X.init = "nndsvd",
    X.L2.ortho = 100.0,
    C1.L1 = 1.0,
    C2.L1 = 0.1,
    epsilon = 1e-6,
    maxit = 5000,
    seed  = 123,
    ...,
    method = c("fiml", "mu"),
    X = NULL,
    C1.restriction = c("union", "none"),
    C1.restriction.threshold = 0.05,
    Phi.restriction = c("full", "diag"),
    C1.L1.path = NULL,
    select = c("BIC", "none")
) {
  cl <- match.call()
  method <- match.arg(method)
  ## `C1.L1` / `C2.L1` penalize the multiplicative updates.  The likelihood path has
  ## its own penalty -- the grid `C1.L1.path`, swept with the support chosen by BIC --
  ## and leaves Theta2 unpenalized, so these two arguments have no meaning under
  ## method = "fiml".  Until 0.9.8 they were accepted there and silently ignored.
  if (method == "fiml" && any(c("C1.L1", "C2.L1") %in% names(cl)))
    warning("`C1.L1` / `C2.L1` apply to method = \"mu\" only and are ignored here. ",
            "Under method = \"fiml\" the L1 penalty on Theta1 is the path `C1.L1.path` ",
            "(swept, with the support selected by BIC) and Theta2 is unpenalized.",
            call. = FALSE)
  ## Options withdrawn in 0.9.8 would otherwise be swallowed by `...` without a word.
  gone <- c(starts = "every penalized fit is warm-started from the unpenalized fit, the only start measurement ever preferred",
            nsplit = "sample splitting was withdrawn: same size as calibration = \"procedure\", lower power",
            calibration = "the calibration level is an argument of nmf.ffb.test(), not of the fit")
  hit <- base::intersect(base::names(gone), base::names(cl))
  if (length(hit) > 0)
    warning("`", paste(hit, collapse = "`, `"), "` was removed in 0.9.8 and is ignored: ",
            gone[[hit[1]]], ".", call. = FALSE)
  ## A RENAMED argument is the more dangerous case: `mask = "none"` reaches `...`, is
  ## dropped, and the fit silently uses the DEFAULT restriction -- the opposite of the
  ## request.  Stop rather than warn.
  renamed <- c(mask = "C1.restriction", cross.threshold = "C1.restriction.threshold",
               phi = "Phi.restriction", lambda1 = "C1.L1.path")
  hit <- base::intersect(base::names(renamed), base::names(cl))
  if (length(hit) > 0)
    stop("`", paste(hit, collapse = "`, `"), "` was renamed in 0.9.8; use `",
         paste(renamed[hit], collapse = "`, `"),
         "`. It would otherwise be dropped silently and the fit would use the default.",
         call. = FALSE)
  if (method == "mu") {
    ## The anchor options exist for the fiml basis only; the legacy estimator
    ## would drop X.anchor without a word, so refuse both.
    if (!is.null(list(...)$X.anchor) || identical(X.init, "spa"))
      stop("nmf.ffb(): X.anchor and X.init = \"spa\" are available with ",
           "method = \"fiml\" only.", call. = FALSE)
    ## Legacy estimator: the body is unchanged (bit-identical) and lives in
    ## .nmf.ffb.mu(); only the `method` field is appended.
    out <- .nmf.ffb.mu(Y1, Y2, rank = rank, X.init = X.init, X.L2.ortho = X.L2.ortho,
                       C1.L1 = C1.L1, C2.L1 = C2.L1, epsilon = epsilon, maxit = maxit,
                       seed = seed, cl = cl, ...)
    out$method <- "mu"
    return(out)
  }
  if (!is.matrix(C1.restriction)) C1.restriction <- match.arg(C1.restriction)
  Phi.restriction <- match.arg(Phi.restriction)
  select <- match.arg(select)
  ## Each penalized fit is warm-started from the unpenalized full-feedback fit and each support refit from
  ## its own penalized fit.  Alternative start sets were measured and never won (NEWS 0.9.8), so there is
  ## no choice left to expose.
  .nmf.ffb.fiml(Y1, Y2, rank = rank, X.init = X.init, X.L2.ortho = X.L2.ortho,
                epsilon = epsilon, maxit = maxit, seed = seed,
                X = X, C1.restriction = C1.restriction, C1.restriction.threshold = C1.restriction.threshold,
                Phi.restriction = Phi.restriction, C1.L1.path = C1.L1.path, select = select, cl = cl, ...)
}

#' Legacy multiplicative-update NMF-FFB estimator (Internal)
#'
#' The body of \code{nmf.ffb()} as it was before \code{method = "fiml"} was
#' added, moved here verbatim so that \code{nmf.ffb(method = "mu")} stays
#' bit-identical to earlier releases.  \code{cl} is the caller's matched call.
#' @keywords internal
#' @noRd
.nmf.ffb.mu <- function(
    Y1, Y2,
    rank = NULL,
    X.init = "nndsvd",
    X.L2.ortho = 100.0,
    C1.L1 = 1.0,
    C2.L1 = 0.1,
    epsilon = 1e-6,
    maxit = 5000,
    seed  = 123,
    cl = NULL,
    ...
) {
  # ------------------------------ checks ------------------------------
  if (!is.matrix(Y1)) Y1 <- as.matrix(Y1)
  if (!is.matrix(Y2)) Y2 <- as.matrix(Y2)

  if (any(!is.finite(Y1)) || any(!is.finite(Y2)))
    stop("Y1 and Y2 must not contain NA/NaN/Inf.")
  if (min(Y1) < 0 || min(Y2) < 0)
    stop("Y1 and Y2 must be non-negative.")
  if (ncol(Y1) != ncol(Y2))
    stop("ncol(Y1) must be equal to ncol(Y2).")

  extra_args <- list(...)
  Q_hidden <- if (!is.null(extra_args$Q)) extra_args$Q else NULL
  Q0 <- if (!is.null(rank)) rank else if (!is.null(Q_hidden)) Q_hidden else nrow(Y2)

  P1 <- nrow(Y1); P2 <- nrow(Y2); N <- ncol(Y1)
  Q  <- Q0
  if (Q < 1) stop("Rank Q must be >= 1.")

  # -------------------------- labels (output) -------------------------
  Y1_labels    <- if (!is.null(rownames(Y1))) rownames(Y1) else paste0("Y1_", 1:P1)
  Y2_labels    <- if (!is.null(rownames(Y2))) rownames(Y2) else paste0("Y2_", 1:P2)
  Basis_labels <- paste0("Factor", 1:Q)

  ## Keep the self-seeding of this fit out of the caller's random stream.
  ## The guard matters: set.seed(NULL) does not mean "do nothing", it
  ## re-initializes the generator from the clock, so an explicit seed = NULL
  ## fit was irreproducible.  With the guard, seed = NULL means "leave the
  ## stream alone" and the fit follows whatever state the caller is in.
  .rng <- .nmfkc.rng.save(seed)
  on.exit(.nmfkc.rng.restore(.rng), add = TRUE)
  if (!is.null(seed)) set.seed(seed)
  .eps <- 1e-10
  .xnorm  <- function(X) sweep(X, 2, pmax(colSums(X), .eps), "/")
  mat1norm <- function(A) max(colSums(abs(A)))

  # ---------- (optional) internal nmfkc warm-start + SC.map baseline --
  ## Resolution rules for `nmfkc.baseline`:
  ##   * not given (default): if X.init is a string method (or NULL),
  ##     run nmfkc(Y1, A = Y2) INTERNALLY using the same X.init,
  ##     X.L2.ortho, epsilon, maxit, seed; use its X as warm-start
  ##     and X * C as M.simple for SC.map.  This is the typical
  ##     workflow described in Satoh (2025) and means
  ##       res <- nmf.ffb(Y1, Y2, rank = Q)
  ##     can run end-to-end without first calling nmfkc().
  ##   * TRUE: same as above (explicit opt-in even when X.init is a
  ##     user-supplied matrix; the matrix is overridden by nmfkc's X).
  ##   * FALSE: opt-out — no internal nmfkc, no warm-start, SC.map = NA
  ##     (equivalent to pre-v0.6.8 behavior).
  ##   * an nmfkc result (list with $X and $C): use as the baseline
  ##     for SC.map and as warm-start for X.init (no internal call).
  baseline_for_scmap <- NULL

  user_baseline    <- extra_args$nmfkc.baseline
  user_M.simple    <- extra_args$M.simple

  baseline_is_obj  <- is.list(user_baseline) &&
                      !is.null(user_baseline$X) && !is.null(user_baseline$C)
  baseline_is_TRUE <- isTRUE(user_baseline)
  baseline_is_FALSE <- isFALSE(user_baseline)
  baseline_is_default <- is.null(user_baseline)

  ## Default auto-run when baseline is unspecified AND X.init is a string
  ## method (NULL is treated as "nndsvd" string).
  auto_nmfkc <- if (baseline_is_FALSE) {
    FALSE
  } else if (baseline_is_TRUE) {
    TRUE
  } else if (baseline_is_obj) {
    FALSE
  } else if (baseline_is_default) {
    is.null(X.init) || is.character(X.init)
  } else {
    FALSE
  }

  if (auto_nmfkc) {
    ## Internal nmfkc call.  Forward only the genuinely shared options
    ## (X.init, X.L2.ortho, epsilon, maxit, seed); the nmf.ffb-specific
    ## C1.L1 / C2.L1 do not apply to the feedforward baseline model.
    nmfkc_xinit <- if (is.null(X.init)) "nndsvd" else X.init
    baseline_for_scmap <- nmfkc(
      Y = Y1, A = Y2, Q = Q,
      X.init = nmfkc_xinit,
      X.L2.ortho = X.L2.ortho,
      epsilon = epsilon,
      maxit = maxit,
      seed = seed,
      verbose = FALSE,
      print.dims = FALSE
    )
    ## Override X.init with the nmfkc-fitted X for nmf.ffb warm-start
    X.init <- baseline_for_scmap$X
  } else if (baseline_is_obj) {
    baseline_for_scmap <- user_baseline
    ## When the user supplies an nmfkc result AND X.init is still a
    ## string / NULL, also use baseline$X as warm-start (rmd workflow).
    if (is.null(X.init) || is.character(X.init)) {
      X.init <- baseline_for_scmap$X
    }
  }

  # ---------------------------- init X,C1,C2 --------------------------
  ## X.init dispatch: delegate string methods to the shared internal
  ## helper .init_X_method() (defined in R/nmfkc.R).  Accepts:
  ##   "nndsvd" (default), "kmeans", "kmeansar", "runif",
  ##   a numeric P1 x Q matrix, or NULL (alias for "nndsvd").
  if (is.null(X.init)) X.init <- "nndsvd"
  if (is.character(X.init)) {
    X <- .init_X_method(X.init, Y1, Q)
  } else {
    X <- as.matrix(X.init)
    if (!all(dim(X) == c(P1, Q))) {
      stop("X.init must have dimension (nrow(Y1) x rank).")
    }
    X[X < 0] <- 0
  }
  X <- .xnorm(X)

  C1 <- matrix(stats::runif(Q * P1, 0.01, 0.1), nrow = Q, ncol = P1)
  C2 <- matrix(stats::runif(Q * P2, 0.001, 0.01), nrow = Q, ncol = P2)

  min_dim <- min(Q, P2)
  for (i in 1:min_dim) {
    C2[i, i] <- stats::runif(1, 0.1, 0.2)
  }

  objfunc      <- numeric(maxit)
  objfunc.penalized <- numeric(maxit)

  # ----------------------------- main loop ----------------------------
  for (it in 1:maxit) {
    M  <- C1 %*% Y1 + C2 %*% Y2

    # 2.1 update X
    Numerator_X       <- tcrossprod(Y1, M)
    Denominator_X_rec <- tcrossprod(X %*% M, M)

    if (X.L2.ortho > 0) {
      XtX <- t(X) %*% X
      XtX_offdiag <- XtX
      diag(XtX_offdiag) <- 0
      Denominator_X_ortho <- X.L2.ortho * X %*% XtX_offdiag
    } else {
      Denominator_X_ortho <- 0
    }

    X <- X * (Numerator_X / (Denominator_X_rec + Denominator_X_ortho + .eps))
    X <- .xnorm(X)

    Xt  <- t(X)
    XtX <- Xt %*% X

    # 2.2 update C1
    XtY1           <- Xt %*% Y1
    Numerator_C1   <- tcrossprod(XtY1, Y1)
    Denominator_C1 <- tcrossprod(XtX %*% M, Y1) + C1.L1 + .eps
    C1 <- C1 * (Numerator_C1 / Denominator_C1)

    # 2.3 update C2
    Numerator_C2   <- tcrossprod(XtY1, Y2)
    Denominator_C2 <- tcrossprod(XtX %*% (C1 %*% Y1 + C2 %*% Y2), Y2) + C2.L1 + .eps
    C2 <- C2 * (Numerator_C2 / Denominator_C2)

    # loss + penalties
    XB <- X %*% (C1 %*% Y1 + C2 %*% Y2)
    loss_rec <- sum((Y1 - XB)^2)
    objfunc[it] <- loss_rec

    if (X.L2.ortho > 0) {
      XtX_off <- XtX
      diag(XtX_off) <- 0
      pen_X_ortho <- 0.5 * X.L2.ortho * sum(XtX_off^2)
    } else {
      pen_X_ortho <- 0
    }
    pen_C1_L1 <- C1.L1 * sum(C1)
    pen_C2_L1 <- C2.L1 * sum(C2)
    objfunc.penalized[it] <- loss_rec + pen_X_ortho + pen_C1_L1 + pen_C2_L1

    if (it >= 10) {
      epsilon_iter <- abs(objfunc[it] - objfunc[it - 1]) / pmax(abs(objfunc[it]), 1)
      if (epsilon_iter <= epsilon) break
    }
  }
  ## Warn when the MU loop exhausts maxit without meeting the
  ## relative-tolerance criterion (matches the nmfkc() convention).
  if (it == maxit && exists("epsilon_iter") && epsilon_iter > abs(epsilon))
    warning(paste0("maximum iterations (", maxit, ") reached..."))

  # ------------------ reorder factors (nmfkc centroid order) ----------
  centroid <- as.numeric((1:nrow(X)) / nrow(X)) %*% X
  index <- order(centroid)
  X  <- X[, index, drop = FALSE]
  C1 <- C1[index, , drop = FALSE]
  C2 <- C2[index, , drop = FALSE]

  # ------------------------------ names -------------------------------
  colnames(X)  <- Basis_labels
  rownames(C1) <- Basis_labels
  rownames(C2) <- Basis_labels
  rownames(X)  <- Y1_labels
  colnames(C1) <- Y1_labels
  colnames(C2) <- Y2_labels

  # -------------------- feedback + stability diagnostics --------------
  XC1  <- X %*% C1
  eigs <- eigen(XC1, only.values = TRUE)$values
  rho  <- max(abs(eigs))
  if (rho >= 1)
    warning("Leontief.inv may be unstable; spectral radius >= 1.")

  XC1_norm1 <- mat1norm(XC1)

  # -------------------- Leontief inverse + equilibrium mapping --------
  I_mat <- diag(nrow(XC1))
  XC2   <- X %*% C2

  Leontief.inv <- tryCatch(
    base::solve(I_mat - XC1),
    error = function(e) {
      warning("Failed to compute Leontief.inv via solve(I - XC1). Returning NA matrices.")
      matrix(NA_real_, nrow = nrow(XC1), ncol = ncol(XC1))
    }
  )
  M.model <- Leontief.inv %*% XC2

  amplification <- mat1norm(M.model) / (mat1norm(XC2) + .eps)
  amplification.bound <- if (XC1_norm1 < 1) 1 / (1 - XC1_norm1) else Inf

  # -------------------- fit indices (equilibrium prediction) ----------
  if (anyNA(Leontief.inv)) {
    Y1_hat <- matrix(NA_real_, nrow = nrow(Y1), ncol = ncol(Y1))
    SC.cov <- NA_real_
    MAE    <- NA_real_
  } else {
    Y1_hat   <- M.model %*% Y2
    S.sample <- Y1 %*% t(Y1)
    S.model  <- Y1_hat %*% t(Y1_hat)
    SC.cov   <- stats::cor(as.numeric(S.sample), as.numeric(S.model))
    MAE      <- mean(abs(Y1 - Y1_hat))
  }

  ## -------------------- input-output structural fidelity (SC.map) -----
  ## SC.map = cor(vec(M.model), vec(M.simple)), where M.simple = X0 * Theta0
  ## is the feedforward baseline mapping (Satoh 2025 §4.SC.map).
  ## `baseline_for_scmap` was set above to either:
  ##   - the result of an internal nmfkc call (default auto path), or
  ##   - the user-supplied `nmfkc.baseline` list, or
  ##   - NULL if the user opted out (nmfkc.baseline = FALSE) or
  ##     supplied X.init as a numeric matrix without nmfkc.baseline.
  ## A user-supplied `M.simple` matrix takes precedence over both.
  SC.map <- NA_real_
  M.simple <- if (!is.null(user_M.simple)) user_M.simple
              else if (!is.null(baseline_for_scmap))
                baseline_for_scmap$X %*% baseline_for_scmap$C
              else NULL
  if (!is.null(M.simple) && !anyNA(M.model)) {
    M.simple <- as.matrix(M.simple)
    if (all(dim(M.simple) == dim(M.model))) {
      SC.map <- tryCatch(
        stats::cor(as.numeric(M.simple), as.numeric(M.model)),
        error = function(e) NA_real_
      )
    } else {
      warning("M.simple has dimension ", paste(dim(M.simple), collapse = "x"),
              " but M.model is ", paste(dim(M.model), collapse = "x"),
              "; SC.map not computed.")
    }
  }

  out <- list(
    ## call / dims match the other fitters so print.nmf() has a header to show.
    call                = cl,
    dims                = sprintf("Y1(%d,%d)~X(%d,%d)[C1(%d,%d)Y1+C2(%d,%d)Y2]",
                                  nrow(Y1), ncol(Y1), nrow(X), ncol(X),
                                  ncol(X), nrow(Y1), ncol(X), nrow(Y2)),
    X                   = X,
    C1                  = C1,
    C2                  = C2,
    XC1                 = XC1,
    XC2                 = XC2,
    XC1.radius          = rho,
    XC1.norm1           = XC1_norm1,
    Leontief.inv        = Leontief.inv,
    M.model             = M.model,
    amplification       = amplification,
    amplification.bound = amplification.bound,
    rank                = Q,
    SC.cov              = SC.cov,
    SC.map              = SC.map,
    mae                 = MAE,
    ## Effective rank of the latent scores B = C1 Y1 + C2 Y2 (Q x N).
    effective.rank      = .effective.rank(C1 %*% Y1 + C2 %*% Y2),
    objfunc             = objfunc[1:it],
    objfunc.penalized        = objfunc.penalized[1:it],
    iter                = it,
    ## Convergence bookkeeping, matching nmfkc / nmfre / nmfae so print.nmf()
    ## and the summaries can say whether the run finished or hit the cap.
    maxit               = maxit,
    epsilon             = epsilon,
    converged           = !(it == maxit && exists("epsilon_iter") &&
                            epsilon_iter > abs(epsilon))
  )
  class(out) <- c("nmf.ffb", "nmf")
  out
}


#' @title Bootstrap inference for NMF-FFB (X held fixed)
#' @description
#' \code{nmf.ffb.inference} performs statistical inference on the structural
#' coefficient matrices \eqn{C_1} (latent feedback, \eqn{\Theta_1}) and
#' \eqn{C_2} (exogenous loading, \eqn{\Theta_2}) from a fitted
#' \code{\link{nmf.ffb}} model, \emph{after} the feed-forward null has been
#' rejected by \code{\link{nmf.ffb.test}}.  The
#' basis \eqn{\hat X} is held fixed throughout, which avoids label switching
#' and gives a clean conditional interpretation: uncertainty of the structural
#' coefficients given the measurement model.  Which resampling scheme is run
#' depends on \code{object$method}.
#'
#' \strong{\code{method = "fiml"}: a parametric bootstrap from the selected
#' model.}  \eqn{B} data sets are drawn from the selected model (with its
#' \eqn{\Theta_1}) and \eqn{(\Theta_1, \Theta_2, \Phi, \psi)} are re-estimated
#' with \eqn{X} and the selected support held fixed.  Confidence intervals are
#' the centred (basic) percentile intervals \eqn{[2\hat\theta -
#' q^*_{1-\alpha/2},\; 2\hat\theta - q^*_{\alpha/2}]}; \code{support_rate} is
#' the share of replicates with \eqn{|\hat\theta^*| >} \code{threshold} (0 for
#' entries outside the selected support), and \code{prob.unsupported} is
#' \eqn{1 - \mathrm{support\_rate}}.
#'
#' These are \emph{post-selection} quantities, conditional on the support that
#' BIC picked, and they can overstate the evidence.  They answer "how large is
#' this entry, given that the procedure kept it", not "is there feedback at
#' all"; the calibrated statement about the presence of feedback is
#' \code{\link{nmf.ffb.test}}, which is the step that belongs before this one.
#' Each replicate seeds itself (\code{seed + B + b}), so the result does not
#' depend on \code{cores}.
#'
#' \strong{\code{method = "mu"} (legacy fits, and objects without a
#' \code{method} field).}  The procedure is a \strong{full pair bootstrap} that holds the basis
#' matrix \eqn{\hat X} from the original fit fixed across all replicates
#' (which avoids label switching and gives a clean conditional
#' interpretation: ``uncertainty of the structural coefficients given the
#' measurement model''):
#' \enumerate{
#'   \item For each replicate \eqn{b = 1, \dots, B}, resample column indices
#'     \eqn{(i_1, \dots, i_N)} with replacement from \eqn{\{1, \dots, N\}}
#'     and form \eqn{Y_1^{(b)} = Y_1[, i]}, \eqn{Y_2^{(b)} = Y_2[, i]}.
#'   \item Re-estimate \eqn{(C_1^{(b)}, C_2^{(b)})} by running the
#'     \code{\link{nmf.ffb}} multiplicative updates \emph{with \eqn{X = \hat X}
#'     held fixed} (no \eqn{X} update; no centroid sort), using the same
#'     \code{C1.L1}, \code{C2.L1} as the original fit.
#'   \item Discard replicates that violate stationarity
#'     (\eqn{\rho(X C_1^{(b)}) \ge 1}) or have an amplification ratio
#'     exceeding the geometric-series bound by more than 1\%.
#' }
#' Because \eqn{C_1, C_2 \ge 0} are non-negative by construction, exact
#' zeros are essentially never observed in the bootstrap distribution.
#' Significance is assessed via a \strong{support rate} at a small display
#' threshold \eqn{\delta} (default \code{0.01}):
#' \deqn{
#'   \mathrm{sup}(c) \;=\; \frac{1}{|\mathrm{valid}|}
#'     \sum_{b \in \mathrm{valid}} \mathbf{1}\!\left( \hat c^{(b)} > \delta \right).
#' }
#' This is a one-sided counterpart of the classical \eqn{p}-value:
#' large \code{support_rate} indicates strong evidence that the entry is
#' meaningfully positive.  Significance markers follow the lavaan
#' convention with the natural correspondence \eqn{p = 1 - \mathrm{sup}}:
#' \code{*} (sup > 0.95), \code{**} (sup > 0.99), \code{***}
#' (sup > 0.999).  Cutoffs use strict greater-than so the rule
#' mirrors the standard R convention for p-values (p < 0.05 / 0.01 /
#' 0.001 → */**/***), translated to support_rate via
#' \eqn{\mathrm{sup} = 1 - p}.
#'
#' @param object A fitted object returned by \code{\link{nmf.ffb}}.  Must
#'   contain \code{X}, \code{C1}, \code{C2}.
#' @param Y1 Endogenous variable matrix (P1 x N).  Must match the data
#'   used in \code{nmf.ffb()}.
#' @param Y2 Exogenous variable matrix (P2 x N).  Same.
#' @param B Number of bootstrap replicates.  Default \code{1000}, the value the
#'   published analysis used; the other inference functions in the package
#'   default their \code{wild.B} to 500, and this one is deliberately left at
#'   1000 so the manuscript results reproduce out of the box.  Note that
#'   \code{***} (sup > 0.999) is only reachable when \emph{every} replicate
#'   clears \code{threshold}, at any \code{B}: raising \code{B} does not add a
#'   finer grade, it makes the same grade stronger evidence (all 1000 rather
#'   than all 500).  Lowering it to 500 roughly halves the running time, since
#'   each replicate is a re-fit.  For \code{method = "fiml"} a replicate is one
#'   FIML fit on the selected support, so \code{B = 1000} is affordable; the
#'   expensive bootstrap is the null one, which \code{\link{nmf.ffb.test}} runs.
#' @param threshold Display threshold \eqn{\delta} for the support rate
#'   \eqn{\Pr_{\mathrm{boot}}(\hat c^{(b)} > \delta)}.  Default
#'   \code{0.01}; entries below this magnitude are treated as effectively
#'   zero in the path diagram.
#' @param boot.level Confidence level for the bootstrap CI.
#'   Default \code{0.95}.
#' @param C1.L1,C2.L1 L1 sparsity penalties used by the original
#'   \code{method = "mu"} fit.  These must match the fit's hyperparameters
#'   for the bootstrap to estimate the correct model.  Defaults
#'   (\code{1.0}, \code{0.1}) match \code{nmf.ffb}'s defaults but you
#'   should pass the actual values used.  Not used for \code{method = "fiml"}.
#' @param seed Base RNG seed for the bootstrap.  For \code{method = "mu"}
#'   each replicate uses \code{seed + b} (resampling) and
#'   \code{seed + 1000 + b} (\eqn{C_1, C_2} initialization); for
#'   \code{method = "fiml"} see Description.  Default \code{123}.
#' @param ... Hidden options.  Shared: \code{cores} (number of parallel
#'   workers, default \code{getOption("mc.cores", 1L)} for \code{"fiml"},
#'   \code{1} for \code{"mu"}; \code{ncores} is accepted as an alias) and
#'   \code{print.trace}.  For \code{method = "fiml"}: \code{factr} and
#'   \code{fiml.maxit} (L-BFGS-B tolerance and cap, default the values
#'   recorded on \code{object}).  The null bootstrap is no longer run here at
#'   all: see \code{\link{nmf.ffb.test}}.  (Removed: \code{boot.null}, which
#'   used to skip it, and the coefficient
#'   intervals).  For \code{method = "mu"}:
#'   \describe{
#'     \item{\code{epsilon}}{Convergence tolerance for the inner fixed-X MU
#'       loop.  Default \code{1e-8} -- deliberately tighter than a plain fit,
#'       because under multiplicative updates an entry heading for the
#'       non-negativity boundary approaches 0 slowly, so a loose tolerance
#'       stops every replicate while such an entry is still above
#'       \code{threshold}.  The support rate is then inflated and significance
#'       over-declared: on a pure-noise design the support of null entries
#'       falls from 0.75 to 0.20 to 0.017 as \code{epsilon} goes 1e-6, 1e-8,
#'       1e-10, while genuinely supported entries stay at 1.  Tighten further
#'       when near-threshold entries matter; loosen only for speed.}
#'     \item{\code{maxit}}{Maximum iterations for the inner MU loop.
#'       Default \code{100000} (the tighter tolerance needs the headroom).}
#'   }
#'
#' @return The input \code{object} with additional bootstrap inference
#'   components:
#' \item{coefficients}{Data frame with rows for every entry of \eqn{C_1}
#'   and \eqn{C_2} and columns \code{Type} ("C1" / "C2"), \code{Basis},
#'   \code{Covariate}, \code{Estimate}, \code{CI_low}, \code{CI_high},
#'   \code{support_rate}, \code{prob.unsupported} (\eqn{= 1 -
#'   \mathrm{support\_rate}}) and \code{sig}.}
#' \item{C1.support.rate, C2.support.rate}{Per-element support rates
#'   (Q x P1 and Q x P2 matrices).}
#' \item{C1.ci.lower, C1.ci.upper, C2.ci.lower, C2.ci.upper}{Per-element
#'   CI bounds (percentile for \code{"mu"}, centred percentile for
#'   \code{"fiml"}).}
#' \item{C1.boot.draws, C2.boot.draws}{Bootstrap distributions: 3D arrays of shape
#'   B x Q x P1 (and B x Q x P2).  Invalid replicates contain \code{NA}.}
#' \item{rho.boot.draws}{Per-replicate spectral radius \eqn{\rho(X C_1^*)}.}
#' \item{AR.boot, iter.boot}{(\code{"mu"} only) per-replicate amplification
#'   ratio and inner-loop iteration count.}
#' #' \item{boot.B, boot.threshold, boot.level}{Inputs recorded for
#'   reproducibility.}
#' \item{boot.n.valid, boot.n.invalid}{Validity counts (for \code{"fiml"}:
#'   of the selected-model bootstrap).}
#' \item{boot.method}{How the replicates were drawn.}
#' \item{identified.approx, identification}{Whether \eqn{(X, [C_1\ C_2])} is
#'   unique up to a permutation of the bases -- the diagnostic of
#'   \code{\link{nmfkc.inference}} with \eqn{A = [Y_1; Y_2]} (see the Details of
#'   \code{\link{nmfkc.ar.latent.inference}}).  The inference is conditional on
#'   \eqn{\hat X}; when \code{TRUE} it concerns the same coefficients whichever
#'   admissible representative is taken.  The exclusion restriction's zeros in
#'   \eqn{C_1} count.  An algebraic check on the point estimate: no re-fits.}
#'
#' The calibrated test of the feed-forward null -- \code{LR.p.boot},
#' \code{prob.select.null}, \code{C1.restriction.change.rate} and the null
#' bootstrap behind them -- is \strong{not} returned here.  It moved to
#' \code{\link{nmf.ffb.test}} in 0.9.8, so that an object carrying coefficient
#' intervals cannot be mistaken for a test of feedback (see NEWS).
#'
#' @section Lifecycle:
#' This function's interface changed at v0.6.8: the legacy 1-step Newton
#' wild bootstrap (with sandwich SE) has been replaced by the full pair
#' bootstrap described above, following the paper revision.  The fields
#' \code{sigma2.used}, \code{C2.se}, \code{C2.se.boot}, \code{C2.p.side}
#' that the previous implementation produced are no longer present.
#'
#' @seealso \code{\link{nmf.ffb}}, \code{\link{nmf.ffb.DOT}}
#' @references
#' Satoh, K. (2025). Applying non-negative matrix factorization with covariates
#'   to structural equation modeling for blind input-output analysis.
#'   arXiv:2512.18250. \url{https://arxiv.org/abs/2512.18250}
#' @export
#' @examples
#' \donttest{
#' Y <- t(iris[, -5])
#' Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ]
#' res  <- nmf.ffb(Y1, Y2, rank = 2)
## The test of the feed-forward null is nmf.ffb.test(); this function gives
#' ## intervals for the entries the test has already licensed reading.
#' inf <- nmf.ffb.inference(res, Y1, Y2, B = 20)   # use B = 1000 in practice
#' head(inf$coefficients)  # estimate, interval and support rate per entry
#'
#' res.mu  <- nmf.ffb(Y1, Y2, rank = 2, method = "mu")
#' res.mu2 <- nmf.ffb.inference(res.mu, Y1, Y2, B = 200)
#' head(res.mu2$coefficients)
#' }
nmf.ffb.inference <- function(object, Y1, Y2,
                               B = 1000L,
                               threshold = 0.01,
                               boot.level = 0.95,
                               C1.L1 = 1.0,
                               C2.L1 = 0.1,
                               seed = 123L,
                               ...) {
  if (is.null(object$X) || is.null(object$C1) || is.null(object$C2))
    stop("object must contain X, C1, and C2 (returned by nmf.ffb).")
  ## Renamed in 0.9.8; through `...` it would be dropped and the default used.
  if ("ci.level" %in% names(match.call()))
    stop("`ci.level` was renamed to `boot.level` in 0.9.8.", call. = FALSE)

  ## Likelihood-based fits take the parametric-bootstrap branch; objects
  ## without a `method` field predate it and are multiplicative-update fits.
  if (identical(object$method, "fiml"))
    return(.ffb.attach.ident(
      .nmf.ffb.inference.fiml(object, Y1, Y2, B = B, threshold = threshold,
                              boot.level = boot.level, seed = seed, what = "ci", ...),
      Y1, Y2, base::list(...)))

  extra_args  <- base::list(...)
  ## Keep our own seeding out of the caller's random stream.
  .rng <- .nmfkc.rng.save(seed)
  on.exit(.nmfkc.rng.restore(.rng), add = TRUE)
  ## Convergence control of the fixed-X re-fits.  Deliberately tighter than a
  ## plain fit: under multiplicative updates a coefficient heading for the
  ## non-negativity boundary approaches 0 slowly, so a loose tolerance stops
  ## every replicate while such an entry is still above the display threshold.
  ## The support rate is then inflated and significance over-declared -- on a
  ## pure-noise design the support of null entries falls 0.75 -> 0.20 -> 0.017
  ## as epsilon goes 1e-6 -> 1e-8 -> 1e-10, while genuinely supported entries
  ## stay at 1.  Tighten further if near-threshold entries matter.
  epsilon     <- if (!is.null(extra_args$epsilon))     extra_args$epsilon     else 1e-8
  maxit       <- if (!is.null(extra_args$maxit))       extra_args$maxit       else 100000L
  ## `cores` is the house name (CONVENTIONS.md 4); `ncores` stays as the alias.
  ncores      <- if (!is.null(extra_args$cores))       extra_args$cores
                 else if (!is.null(extra_args$ncores)) extra_args$ncores else 1L
  print.trace <- if (!is.null(extra_args$print.trace)) extra_args$print.trace else FALSE

  Y1 <- base::as.matrix(Y1)
  Y2 <- base::as.matrix(Y2)
  if (ncol(Y1) != ncol(Y2)) stop("Y1 and Y2 must have the same number of columns.")

  X.hat  <- object$X
  C1.hat <- object$C1
  C2.hat <- object$C2

  Q  <- ncol(X.hat)
  P1 <- nrow(Y1)
  P2 <- nrow(Y2)
  N  <- ncol(Y1)
  if (nrow(X.hat) != P1)
    stop("nrow(X) must equal nrow(Y1); did Y1 change shape?")
  if (Q < 1L) stop("Inferred rank Q < 1.")

  ## ----------------------------------------------------------------
  ## one_boot: a single bootstrap replicate.
  ## Closure captures Y1, Y2, X.hat, Xt, XtX, C1.L1, C2.L1, Q, P1, P2,
  ## N, epsilon, maxit, seed.  Returns a list with valid flag plus
  ## (C1, C2, rho, AR, iter) when valid, or NA placeholders otherwise.
  ## X.hat is fixed by design, so Xt / XtX are replicate-invariant and
  ## computed once here rather than inside every replicate.
  ## ----------------------------------------------------------------
  Xt  <- t(X.hat)
  XtX <- Xt %*% X.hat
  one_boot <- function(b) {
    ## Sample column indices with replacement (pair bootstrap)
    set.seed(seed + b)
    idx  <- sample.int(N, size = N, replace = TRUE)
    Y1_b <- Y1[, idx, drop = FALSE]
    Y2_b <- Y2[, idx, drop = FALSE]

    ## Independent random init for C1, C2 (X is held at X.hat)
    set.seed(seed + 1000L + b)
    C1 <- matrix(stats::runif(Q * P1, 0.01, 0.1),  nrow = Q)
    C2 <- matrix(stats::runif(Q * P2, 0.001, 0.01), nrow = Q, ncol = P2)
    min_dim <- min(Q, P2)
    if (min_dim >= 1L) {
      for (i in 1:min_dim) C2[i, i] <- stats::runif(1, 0.1, 0.2)
    }

    .eps_local <- 1e-10
    XtY1_b  <- Xt %*% Y1_b              # Q x N (re-used via crossprods below)
    XtY1Y1t <- tcrossprod(XtY1_b, Y1_b) # Q x P1   = Xt %*% Y1_b %*% t(Y1_b)
    XtY1Y2t <- tcrossprod(XtY1_b, Y2_b) # Q x P2
    prev_loss <- Inf
    iter_used <- 0L

    ## Fixed-X MU loop on (C1, C2).  M_b = C1 Y1_b + C2 Y2_b is computed
    ## before the loop and carried across iterations: the value formed after
    ## the C2 update (used for the loss) is exactly the loop-top value of the
    ## next iteration, so each iteration computes it twice instead of three
    ## times with identical operands.
    M_b <- C1 %*% Y1_b + C2 %*% Y2_b     # Q x N latent representation
    for (it in seq_len(maxit)) {
      ## C1 update
      Den_C1 <- tcrossprod(XtX %*% M_b, Y1_b) + C1.L1 + .eps_local
      C1 <- C1 * (XtY1Y1t / Den_C1)

      ## C2 update (uses updated C1 via fresh M_b)
      M_b <- C1 %*% Y1_b + C2 %*% Y2_b
      Den_C2 <- tcrossprod(XtX %*% M_b, Y2_b) + C2.L1 + .eps_local
      C2 <- C2 * (XtY1Y2t / Den_C2)

      ## Convergence check (relative change in reconstruction loss)
      M_b <- C1 %*% Y1_b + C2 %*% Y2_b   # new C1, new C2; carried to next iter
      XB <- X.hat %*% M_b
      loss <- sum((Y1_b - XB)^2)
      if (it >= 10L) {
        rel <- abs(loss - prev_loss) / max(abs(loss), 1)
        if (rel <= epsilon) { iter_used <- it; break }
      }
      prev_loss <- loss
      iter_used <- it
    }

    ## ---- Validity checks (mirrors the paper's bootstrap helper) ----
    XC1_b <- X.hat %*% C1
    rho <- tryCatch(
      max(abs(eigen(XC1_b, only.values = TRUE, symmetric = FALSE)$values)),
      error = function(e) NA_real_
    )
    norm1_XC1 <- max(colSums(abs(XC1_b)))
    AR_bound <- if (is.finite(norm1_XC1) && norm1_XC1 < 1)
                  1 / (1 - norm1_XC1)
                else Inf

    if (!is.finite(rho) || rho >= 1) {
      return(list(valid = FALSE, C1 = NULL, C2 = NULL,
                  rho = rho, AR = NA_real_, iter = iter_used))
    }
    ## XC1_b = X (P1 x Q) %*% C1 (Q x P1) is P1 x P1, so the identity
    ## here must be P1 x P1 too.  (Earlier draft used diag(Q) which is
    ## conformable only when Q == P1, otherwise solve() throws an error
    ## that the tryCatch swallowed -- causing every replicate to be
    ## marked invalid with AR = NA.)
    I_mat   <- diag(nrow(XC1_b))
    Linv    <- tryCatch(solve(I_mat - XC1_b), error = function(e) NULL)
    if (is.null(Linv)) {
      return(list(valid = FALSE, C1 = NULL, C2 = NULL,
                  rho = rho, AR = NA_real_, iter = iter_used))
    }
    XC2_b   <- X.hat %*% C2
    M.model <- Linv %*% XC2_b
    AR <- max(colSums(abs(M.model))) / (max(colSums(abs(XC2_b))) + .eps_local)
    if (!is.finite(AR) || AR > AR_bound * 1.01) {
      return(list(valid = FALSE, C1 = NULL, C2 = NULL,
                  rho = rho, AR = AR, iter = iter_used))
    }

    list(valid = TRUE, C1 = C1, C2 = C2,
         rho = rho, AR = AR, iter = iter_used)
  }

  ## ----------------------------------------------------------------
  ## Run the B replicates (parallel on Linux/macOS via mclapply,
  ## PSOCK cluster on Windows; serial when ncores == 1).
  ## ----------------------------------------------------------------
  if (print.trace)
    base::message(sprintf("  Bootstrap: B=%d, ncores=%d, threshold=%.3g, boot.level=%.2f",
                          B, ncores, threshold, boot.level))

  if (ncores > 1L) {
    if (.Platform$OS.type == "windows") {
      cl <- parallel::makeCluster(ncores)
      on.exit(parallel::stopCluster(cl), add = TRUE)
      parallel::clusterExport(cl,
        varlist = c("Y1", "Y2", "X.hat", "Xt", "XtX", "C1.L1", "C2.L1",
                    "Q", "P1", "P2", "N", "epsilon", "maxit", "seed"),
        envir = environment())
      res_list <- parallel::parLapply(cl, seq_len(B), one_boot)
    } else {
      res_list <- parallel::mclapply(seq_len(B), one_boot, mc.cores = ncores)
    }
  } else {
    res_list <- lapply(seq_len(B), one_boot)
  }

  ## ----------------------------------------------------------------
  ## Aggregate replicates into 3D arrays (B x Q x P1) and (B x Q x P2)
  ## ----------------------------------------------------------------
  C1.boot.draws <- array(NA_real_, dim = c(B, Q, P1))
  C2.boot.draws <- array(NA_real_, dim = c(B, Q, P2))
  rho.vec  <- rep(NA_real_, B)
  AR.vec   <- rep(NA_real_, B)
  iter.vec <- rep(NA_integer_, B)
  valid.vec <- logical(B)

  for (b in seq_len(B)) {
    r <- res_list[[b]]
    if (is.null(r) || inherits(r, "try-error")) next
    valid.vec[b] <- isTRUE(r$valid)
    if (!is.null(r$rho))  rho.vec[b]  <- r$rho
    if (!is.null(r$AR))   AR.vec[b]   <- r$AR
    if (!is.null(r$iter)) iter.vec[b] <- as.integer(r$iter)
    if (valid.vec[b]) {
      C1.boot.draws[b, , ] <- r$C1
      C2.boot.draws[b, , ] <- r$C2
    }
  }
  n.valid <- sum(valid.vec)

  if (n.valid < 10L) {
    warning(sprintf("Only %d / %d bootstrap replicates were valid; CIs / support rates may be unreliable. Consider checking convergence (raise maxit, lower epsilon) or the stationarity of the original fit (rho < 1).",
                    n.valid, B))
  }

  ## ----------------------------------------------------------------
  ## Per-element summary: support rate at threshold, percentile CI
  ## ----------------------------------------------------------------
  alpha <- 1 - boot.level

  apply_finite <- function(arr, FUN) {
    apply(arr, c(2, 3), function(v) {
      v <- v[is.finite(v)]
      if (!length(v)) NA_real_ else FUN(v)
    })
  }

  C1.support  <- apply_finite(C1.boot.draws, function(v) mean(v > threshold))
  C2.support  <- apply_finite(C2.boot.draws, function(v) mean(v > threshold))
  C1.ci.lower <- apply_finite(C1.boot.draws, function(v) stats::quantile(v, alpha / 2,    names = FALSE))
  C1.ci.upper <- apply_finite(C1.boot.draws, function(v) stats::quantile(v, 1 - alpha / 2, names = FALSE))
  C2.ci.lower <- apply_finite(C2.boot.draws, function(v) stats::quantile(v, alpha / 2,    names = FALSE))
  C2.ci.upper <- apply_finite(C2.boot.draws, function(v) stats::quantile(v, 1 - alpha / 2, names = FALSE))

  ## Significance markers from support rate (one-sided; lavaan / stats
  ## convention).  Cutoffs use strict greater-than so the rule mirrors
  ## the R convention for p-values (p < 0.05 / 0.01 / 0.001 → */**/***),
  ## translated to support_rate via support_rate = 1 - p.  This matches
  ## the verbal description "support rates exceeding 0.95, 0.99, 0.999"
  ## in Satoh (2025), arXiv:2512.18250.
  sig.from.support <- function(s) {
    ifelse(!is.finite(s), " ",
      ifelse(s > 0.999, "***",
        ifelse(s > 0.99,  "**",
          ifelse(s > 0.95, "*", " "))))
  }

  ## ----------------------------------------------------------------
  ## Build coefficients table (rows for both C1 and C2)
  ##
  ## `1 - support_rate` is the fraction of bootstrap replicates in which the
  ## coefficient was NOT supported.  That is a tail proportion, not a p-value,
  ## and CONVENTIONS.md section 6 forbids naming such a quantity `p.*` for
  ## exactly this reason -- calling it p_value invites reading 0.03 as
  ## significance at the 5% level, which it is not.  The column is therefore
  ## now `prob.unsupported`, with `p_value` kept alongside it as a
  ## back-compatible duplicate: nmf.ffb.DOT and any user script that filters or
  ## stars on p_value keep working unchanged.
  ## ----------------------------------------------------------------
  Q_lab  <- if (!is.null(rownames(C1.hat))) rownames(C1.hat) else paste0("Factor", 1:Q)
  Y1_lab <- if (!is.null(colnames(C1.hat))) colnames(C1.hat) else paste0("Y1_", 1:P1)
  Y2_lab <- if (!is.null(colnames(C2.hat))) colnames(C2.hat) else paste0("Y2_", 1:P2)

  build_block <- function(type_label, Mhat, sup, lo, hi, basislabs, varlabs) {
    Q_local <- nrow(Mhat); P_local <- ncol(Mhat)
    s_vec <- as.vector(sup)
    data.frame(
      Type         = type_label,
      Basis        = rep(basislabs, times = P_local),
      Covariate    = rep(varlabs,   each  = Q_local),
      Estimate     = as.vector(Mhat),
      CI_low       = as.vector(lo),
      CI_high      = as.vector(hi),
      support_rate = s_vec,
      ## the convention-compliant name ...
      prob.unsupported = ifelse(is.finite(s_vec), 1 - s_vec, NA_real_),
      ## ... and the historical one, same numbers, kept for back-compatibility
      sig          = sig.from.support(s_vec),
      stringsAsFactors = FALSE
    )
  }

  C1_block <- build_block("C1", C1.hat, C1.support, C1.ci.lower, C1.ci.upper, Q_lab, Y1_lab)
  C2_block <- build_block("C2", C2.hat, C2.support, C2.ci.lower, C2.ci.upper, Q_lab, Y2_lab)
  coefficients <- rbind(C1_block, C2_block)
  rownames(coefficients) <- NULL

  if (print.trace)
    base::message(sprintf("  Bootstrap done: %d / %d valid replicates.", n.valid, B))

  ## ----------------------------------------------------------------
  ## Append to object and return
  ## ----------------------------------------------------------------
  object$boot.B          <- B
  object$boot.threshold  <- threshold
  object$boot.level           <- boot.level
  object$boot.n.valid    <- n.valid
  object$boot.n.invalid  <- B - n.valid
  object$rho.boot.draws             <- rho.vec
  object$AR.boot              <- AR.vec
  object$iter.boot            <- iter.vec
  object$C1.boot.draws             <- C1.boot.draws
  object$C2.boot.draws             <- C2.boot.draws
  object$C1.support.rate      <- C1.support
  object$C2.support.rate      <- C2.support
  object$C1.ci.lower          <- C1.ci.lower
  object$C1.ci.upper          <- C1.ci.upper
  object$C2.ci.lower          <- C2.ci.lower
  object$C2.ci.upper          <- C2.ci.upper
  object$coefficients         <- coefficients
  object <- .ffb.attach.ident(object, Y1, Y2, extra_args)

  ## Final class vector for a typical input:
  ##   c("nmf.ffb.inference", "nmf.inference", "nmf.ffb", "nmf")
  if (!inherits(object, "nmf.ffb.inference"))
    ## "nmf.inference" too, so print() reaches print.nmf.inference -> the
    ## coefficient table.  Without it the class chain ended at "nmf" and
    ## printing an inference result showed only the fit header.
    class(object) <- c("nmf.ffb.inference", "nmf.inference", class(object))
  object
}


#' @title Cross-Validation for NMF-FFB
#' @description
#' Performs K-fold cross-validation to evaluate the equilibrium mapping of
#' the NMF-FFB model.
#'
#' For each fold, \code{nmf.ffb} is fitted on the training samples,
#' yielding an equilibrium mapping \eqn{\hat Y_1 = M_{\mathrm{model}} Y_2}.
#' The held-out endogenous variables \eqn{Y_1} are then predicted from \eqn{Y_2}
#' using this mapping, and the mean absolute error (MAE) over all entries in the
#' test block is computed. The returned value is the average MAE across folds.
#'
#' This implements the hyperparameter selection strategy described in the paper:
#' hyperparameters are chosen by predictive cross-validation rather than direct
#' inspection of the internal structural matrices.
#'
#' \strong{\code{method = "fiml"}.}  When \code{method = "fiml"} is passed
#' (through \code{...}), the column-wise scheme above is not run: the
#' equilibrium prediction \eqn{M_{\mathrm{model}} Y_2} is a reduced-form
#' quantity, and the reduced form is the same with and without feedback, so
#' predicting held-out columns cannot select \eqn{\Theta_1} (that is what the
#' BIC path and the bootstrap-calibrated LR test in \code{\link{nmf.ffb}} /
#' \code{\link{nmf.ffb.inference}} are for).  The only tunable quantity is the
#' stage-1 rank, and the call is delegated to the element-wise CV of the
#' feed-forward fit, \code{\link{nmfkc.ecv}(Y = Y1, A = Y2, rank = rank,
#' ...)}, whose object is returned (\code{$objfunc} holds the CV error per
#' rank; \code{rank} may be a vector, default \code{1:min(3, P1)}).
#' \code{C1.L1} and \code{C2.L1} are ignored in that case.
#'
#' @param Y1 A non-negative numeric matrix of endogenous variables with
#'   \strong{rows = variables (P1), columns = samples (N)}.
#' @param Y2 A non-negative numeric matrix of exogenous variables with
#'   \strong{rows = variables (P2), columns = samples (N)}.
#'   Must satisfy \code{ncol(Y1) == ncol(Y2)}.
#' @param rank Integer; rank (number of latent factors) passed to \code{nmf.ffb}.
#'   If \code{NULL}, \code{nmf.ffb} decides the effective rank (via \code{...} or \code{nrow(Y2)}).
#' @param X.init Initialization strategy for \code{X}, forwarded to
#'   \code{\link{nmf.ffb}}.  One of \code{"nndsvd"} (default),
#'   \code{"kmeans"}, \code{"kmeansar"}, \code{"runif"}, a numeric
#'   \eqn{P_1 \times Q} matrix, or \code{NULL} (alias for
#'   \code{"nndsvd"}).  See \code{\link{nmf.ffb}} for details.
#' @param X.L2.ortho L2 orthogonality penalty for \code{X}.
#' @param C1.L1 L1 sparsity penalty for \code{C1} (\eqn{\Theta_1}).
#' @param C2.L1 L1 sparsity penalty for \code{C2} (\eqn{\Theta_2}).
#' @param epsilon Convergence threshold for \code{nmf.ffb}.
#' @param maxit Maximum number of iterations for \code{nmf.ffb}.
#' @param ... Additional arguments passed to \code{nmf.ffb} (except for
#'   \code{rank}, \code{seed}, \code{div}, \code{shuffle}, which are handled here).
#'   Also accepts: \code{nfolds} (number of folds, default 5; \code{div} also accepted),
#'   \code{seed} (master random seed, default \code{NULL}),
#'   \code{shuffle} (logical, default \code{TRUE}),
#'   \code{cores} (integer; number of parallel workers for the fold loop,
#'   default \code{getOption("mc.cores", 1L)}). Fold partitions and per-fold
#'   seeds are drawn before the loop and each \code{nmf.ffb} fit self-seeds,
#'   so results are identical to the sequential run for any \code{cores}.
#'
#' @return A numeric scalar: mean MAE across CV folds (the fits use
#'   \code{method = "mu"} unless \code{method} is given).  With
#'   \code{method = "fiml"}, the \code{\link{nmfkc.ecv}} object of the
#'   stage-1 rank sweep.
#'
#' @examples
#' Y <- t(iris[, -5])
#' Y1 <- Y[1:2, ]
#' Y2 <- Y[3:4, ]
#' mae <- nmf.ffb.cv(Y1, Y2, rank = 2, maxit = 500, nfolds = 3)
#' mae
#' \donttest{
#' # rank selection for the likelihood-based estimator
#' ecv <- nmf.ffb.cv(Y1, Y2, rank = 1:2, method = "fiml", nfolds = 3)
#' ecv$objfunc
#' }
#'
#' @seealso \code{\link{nmf.ffb}}
#' @export
nmf.ffb.cv <- function(
    Y1, Y2,
    rank = NULL,
    X.init = "nndsvd",
    X.L2.ortho = 100.0,
    C1.L1 = 1.0,        # L1 sparsity for C1 (Theta1)
    C2.L1 = 0.1,        # L1 sparsity for C2 (Theta2)
    epsilon = 1e-6,     # Convergence tolerance passed to nmf.ffb
    maxit = 5000,
    ...
){
  extra_cv <- base::list(...)
  ## method = "fiml": column-wise CV of the equilibrium mapping cannot select
  ## Theta1 (the reduced form M Y2 is the same with and without feedback), so
  ## the only tunable quantity is the stage-1 rank; delegate to the
  ## element-wise CV of the feed-forward fit.
  if (identical(extra_cv[["method", exact = TRUE]], "fiml")) {
    if (is.null(rank)) rank <- seq_len(min(3L, nrow(as.matrix(Y1))))
    ecv_args <- extra_cv
    ecv_args$method <- NULL
    ecv_args$C1.L1 <- NULL; ecv_args$C2.L1 <- NULL
    for (nm in c("X", "C1.restriction", "C1.restriction.threshold", "Phi.restriction",
                 "C1.L1.path", "select",
                 "fiml.maxit", "factr", "nmfkc.baseline", "M.simple"))
      ecv_args[[nm]] <- NULL
    if (is.null(ecv_args$X.init)) ecv_args$X.init <- X.init
    return(do.call("nmfkc.ecv", c(list(Y = as.matrix(Y1), A = as.matrix(Y2), rank = rank,
                                       X.L2.ortho = X.L2.ortho, epsilon = epsilon,
                                       maxit = maxit), ecv_args)))
  }
  nfolds  <- if (!is.null(extra_cv$nfolds))  extra_cv$nfolds  else if (!is.null(extra_cv$div)) extra_cv$div else 5
  ## CONVENTIONS.md 2: default seeds are concrete (123), so an analysis is
  ## reproducible without the user having to think about it.  This was the one
  ## CV in the package defaulting to NULL, i.e. the one whose fold assignment
  ## changed between runs.  Pass seed = NULL explicitly for the old behaviour.
  seed    <- if ("seed" %in% names(extra_cv)) extra_cv$seed else 123L
  ## Keep our own seeding out of the caller's random stream (no-op when
  ## seed = NULL, which here means "use the stream as it is").
  .rng <- .nmfkc.rng.save(seed)
  on.exit(.nmfkc.rng.restore(.rng), add = TRUE)
  shuffle <- if (!is.null(extra_cv$shuffle)) extra_cv$shuffle else TRUE
  cores   <- if (!is.null(extra_cv$cores))   extra_cv$cores   else getOption("mc.cores", 1L)
  div <- nfolds
  # ------------------------------------------------------------------
  # 1. Basic input checks
  #
  # NMF-FFB requires non-negative matrices. We also require that Y1 and Y2
  # share the same number of samples (columns) to allow paired CV splits.
  # ------------------------------------------------------------------
  if (!is.matrix(Y1)) Y1 <- as.matrix(Y1)
  if (!is.matrix(Y2)) Y2 <- as.matrix(Y2)
  if (any(!is.finite(Y1)) || any(!is.finite(Y2))) stop("Y1 and Y2 must not contain NA/NaN/Inf.")
  div <- as.integer(div)
  if (min(Y1) < 0 || min(Y2) < 0) stop("Y1 and Y2 must be non-negative.")
  if (ncol(Y1) != ncol(Y2)) {
    stop("Y1 and Y2 must have the same number of columns (samples).")
  }

  P1 <- nrow(Y1)
  P2 <- nrow(Y2)
  N  <- ncol(Y1)

  if (div < 2L) {
    stop("div (number of CV folds) must be >= 2.")
  }
  if (div > N) {
    stop("div (number of CV folds) must be <= number of samples.")
  }

  # ------------------------------------------------------------------
  # 2. Handle extra arguments for nmf.ffb
  #
  # We collect additional arguments in 'extra_args' but explicitly remove
  # those that are managed at the CV level:
  #   - div, shuffle : used only here, not passed to nmf.ffb.
  #   - rank        : passed explicitly from nmf.ffb.cv.
  #   - seed        : fold-specific seeds are generated here.
  # ------------------------------------------------------------------
  extra_args <- list(...)

  extra_args$div     <- NULL
  extra_args$nfolds  <- NULL
  extra_args$shuffle <- NULL
  extra_args$rank    <- NULL
  extra_args$seed    <- NULL
  extra_args$cores   <- NULL
  ## The column-wise scheme tunes C1.L1 / C2.L1, which only the
  ## multiplicative-update estimator uses, so the fold fits stay on
  ## method = "mu" unless the caller says otherwise (a "fiml" request was
  ## delegated to nmfkc.ecv above).  This keeps the CV scores identical to
  ## those before "fiml" became nmf.ffb()'s default.
  if (is.null(extra_args[["method", exact = TRUE]])) extra_args$method <- "mu"

  # ------------------------------------------------------------------
  # 3. Set RNG for CV partition and per-fold seeds
  #
  # If a master 'seed' is given:
  #   - it is used to define the CV partition (sample permutation),
  #   - independent seeds for each nmf.ffb run are drawn.
  # If 'seed' is NULL:
  #   - CV partition uses the current RNG state,
  #   - nmf.ffb runs use whatever the global RNG state is at call time.
  # ------------------------------------------------------------------
  if (!is.null(seed)) {
    set.seed(seed)
  }

  # Sample indices for fold division
  if (shuffle) {
    perm_index <- sample.int(N)
  } else {
    perm_index <- seq_len(N)
  }

  # Per-fold seeds for nmf.ffb (optional, only if master seed specified)
  if (!is.null(seed)) {
    seeds_fold <- sample.int(.Machine$integer.max, div)
  } else {
    seeds_fold <- rep(NA_integer_, div)
  }

  # ------------------------------------------------------------------
  # 4. Create CV folds
  #
  # We assign approximately N/div samples to each fold, distributing
  # any remainder one-by-one to the earliest folds.
  # 'block[i] = k' means sample i belongs to fold k (as test set).
  # ------------------------------------------------------------------
  remainder <- N %% div
  division  <- N %/% div
  block     <- integer(N)

  processed_count <- 0L
  for (i in 1:(div - 1L)) {
    plus       <- ifelse(i <= remainder, 1L, 0L)
    chunk_size <- division + plus
    idx_range  <- (processed_count + 1L):(processed_count + chunk_size)
    target_idx <- perm_index[idx_range]
    block[target_idx] <- i
    processed_count   <- processed_count + chunk_size
  }
  # Last fold gets all remaining samples
  target_idx <- perm_index[(processed_count + 1L):N]
  block[target_idx] <- div

  # Per-fold CV loss (MAE) will be stored here
  objfunc.block <- numeric(div)

  # ------------------------------------------------------------------
  # 5. Cross-validation loop
  #
  # For each fold j:
  #   - train on all samples not in fold j,
  #   - test on samples in fold j,
  #   - fit nmf.ffb on training data,
  #   - compute MAE on test block from equilibrium mapping M.model.
  # ------------------------------------------------------------------
  # Per-fold task. The partition 'block' and per-fold seeds 'seeds_fold' are
  # fixed above, and the only RNG consumer here is nmf.ffb() (self-seeds), so
  # folds are independent and deterministic; .nmfkc.parlapply returns results
  # in input (fold) order, so cores > 1 is bit-identical to the sequential run.
  fold_fun <- function(j) {
    # Train / Test split
    train_idx <- block != j
    test_idx  <- block == j

    Y1_train <- Y1[, train_idx, drop = FALSE]
    Y1_test  <- Y1[, test_idx,  drop = FALSE]
    Y2_train <- Y2[, train_idx, drop = FALSE]
    Y2_test  <- Y2[, test_idx,  drop = FALSE]

    # Fold-specific seed for nmf.ffb
    seed_j <- if (!is.null(seed)) seeds_fold[j] else NULL

    # Assemble arguments for nmf.ffb
    ffb_args <- c(
      extra_args,   # User-specified additional arguments (e.g., Q)
      list(
        Y1         = Y1_train,
        Y2         = Y2_train,
        rank       = rank,
        X.init     = X.init,
        X.L2.ortho = X.L2.ortho,
        C1.L1      = C1.L1,
        C2.L1      = C2.L1,
        epsilon    = epsilon,
        maxit      = maxit
      )
    )
    # Attach seed only when it is defined
    if (!is.null(seed_j)) {
      ffb_args$seed <- seed_j
    }

    # Call nmf.ffb on the training data (suppress messages for cleaner CV output)
    res_j <- suppressMessages(do.call("nmf.ffb", ffb_args))

    # If mapping is not usable, penalize this fold (do not crash CV)
    if (is.null(res_j$M.model) || any(!is.finite(res_j$M.model))) {
      return(Inf)
    }
    Pre_test <- res_j$M.model %*% Y2_test
    if (any(!is.finite(Pre_test))) {
      return(Inf)
    }
    mean(abs(Y1_test - Pre_test))
  }
  objfunc.block <- base::unlist(
    .nmfkc.parlapply(base::as.list(1:div), fold_fun,
                     cores = cores, envir = environment()),
    use.names = FALSE)

  # ------------------------------------------------------------------
  # 6. Aggregate CV score
  #
  # The overall CV criterion is the average MAE across all folds.
  # This is typically minimized over hyperparameter grids
  # (e.g., rank, X.L2.ortho, C1.L1, C2.L1) when tuning NMF-FFB.
  # ------------------------------------------------------------------
  objfunc <- mean(objfunc.block)
  return(objfunc)
}



#' @title Structural Diagnostics of a Fitted NMF-FFB Feedback Matrix
#'
#' @description
#' Two quantities that describe \emph{what} the selected feedback is, as opposed
#' to whether it is there (which is the job of the calibrated likelihood-ratio
#' test in \code{\link{nmf.ffb.inference}}).
#'
#' \strong{Cycle structure.} The equilibrium operator is \eqn{A = X\Theta_1} on
#' the indicators and \eqn{A_Q = \Theta_1X} on the factors, with the same
#' spectral radius \eqn{\rho}.  The strongly connected components of the
#' directed graph of \eqn{A_Q} are the cycles.  \eqn{\rho = 0} with a non-empty
#' support means \eqn{A} is nilpotent: the selected feedback is a cascade of
#' directed cross-factor paths without a closed loop.  Since \eqn{A \ge 0},
#' Perron--Frobenius gives a non-negative leading eigenvector, and its
#' normalized entries say how much of the loop each factor carries.
#'
#' \strong{Factor-level alignment.} Feedback of the form
#' \eqn{\Theta_1 = aX^\top\Psi^{-1}} enters each factor through the factor
#' scores of the factors and is observationally equivalent, under the Gaussian
#' working model, to a change of \eqn{\Theta_2} and of the factor covariance
#' \eqn{\Phi}; only the departure from that family is identified, through the
#' uniquenesses \eqn{\Psi}.  On the free set \eqn{F} of the exclusion restriction the
#' family is a subspace of dimension at most \eqn{Q(Q-1)}, and
#' \code{omega} is the share of the squared norm of \eqn{\Theta_1} that lies in
#' it.  \code{omega} has no natural zero: for an isotropic direction its
#' expectation is \code{omega0 = dim / n.free}, so it should be read against
#' \code{omega0} and against a simulation in which the true feedback is
#' indicator-level (there it averages about \code{1.3 * omega0}).
#'
#' @param object A fitted \code{nmf.ffb} object from
#'   \code{\link{nmf.ffb}} with \code{method = "fiml"}.
#' @param which Which feedback matrix to describe: \code{"selected"} (default,
#'   the BIC-selected support) or \code{"full"} (the unpenalized fit).
#' @return A list with \code{cycles} (\code{A.factor}, \code{rho},
#'   \code{cycle}, \code{n.cycles}, \code{perron}) and \code{omega}
#'   (\code{omega}, \code{omega0}, \code{dim}, \code{n.free}).
#' @seealso \code{\link{nmf.ffb}}, \code{\link{nmf.ffb.inference}}
#' @examples
#' \donttest{
#' set.seed(1)
#' Y2 <- matrix(runif(2 * 200), 2, 200)
#' X <- matrix(0.02, 6, 2); X[1:3, 1] <- 1; X[4:6, 2] <- 1
#' X <- sweep(X, 2, colSums(X), "/")
#' T1 <- matrix(0, 2, 6); T1[1, 5] <- 0.6; T1[2, 2] <- 0.5
#' L <- solve(diag(6) - X %*% T1)
#' Y1 <- pmax(L %*% (X %*% (matrix(c(1, 0.5, 0.3, 1), 2, 2) %*% Y2) +
#'                   matrix(rnorm(6 * 200, 0, 0.05), 6, 200)), 0)
#' fit <- nmf.ffb(Y1, Y2, Q = 2)
#' d <- nmf.ffb.diagnostics(fit)
#' d$cycles$rho          # spectral radius
#' d$cycles$cycle        # which factors lie on a cycle
#' c(d$omega$omega, d$omega$omega0)
#' }
#' @export
nmf.ffb.diagnostics <- function(object, which = c("selected", "full")) {
  which <- base::match.arg(which)
  if (!base::inherits(object, "nmf.ffb"))
    base::stop("nmf.ffb.diagnostics() expects an object from nmf.ffb().")
  if (base::is.null(object$C1.free) || base::is.null(object$psi))
    base::stop("this fit carries no exclusion restriction / uniquenesses: nmf.ffb.diagnostics() needs method = \"fiml\".")
  T1 <- if (which == "selected") object$C1 else object$full$C1
  psi <- if (which == "selected") object$psi else object$full$psi
  if (base::is.null(T1)) base::stop("no ", which, " feedback matrix on this object.")
  ## The best few distinct supports, so that the report says how thin the margin is.  BIC selects one
  ## support; adding a free entry costs log(N), and a difference in BIC is twice the log of a Bayes factor,
  ## so a difference below about 2 is not evidence for one support over the other (Kass and Raftery 1995).
  ## When the best supports form a chain under inclusion their intersection is the part the criterion does
  ## not put in doubt; when they do not, no single support should be read off.
  top <- NULL; core <- NULL; envelope <- NULL; nested <- NA
  if (!base::is.null(object$candidates) && !base::is.null(object$supports)) {
    cd <- object$candidates[base::order(object$candidates$BIC, object$candidates$nnz), , drop = FALSE]
    k <- base::min(3L, base::nrow(cd))
    top <- base::data.frame(order = base::seq_len(k), nnz = cd$nnz[base::seq_len(k)],
                            rho = cd$rho[base::seq_len(k)], BIC = cd$BIC[base::seq_len(k)],
                            dBIC = cd$BIC[base::seq_len(k)] - cd$BIC[1],
                            stringsAsFactors = FALSE)
    tops <- base::lapply(base::seq_len(k), function(i) object$supports[[cd$support_id[i]]])
    core <- base::Reduce(`&`, tops) * 1
    envelope <- base::Reduce(`|`, tops) * 1
    nested <- base::all(base::outer(base::seq_len(k), base::seq_len(k), base::Vectorize(function(a, b)
      base::all(!tops[[a]] | tops[[b]]) || base::all(!tops[[b]] | tops[[a]]))))
  }
  base::list(cycles = .ffb.fiml.cycles(object$X, T1),
             omega = .ffb.fiml.omega(object$X, T1, psi, object$C1.free),
             top = top, core = core, envelope = envelope, nested = nested,
             log.N = if (!base::is.null(object$N)) base::log(object$N) else NA_real_,
             which = which)
}

#' @title Heuristic Variable Splitting for NMF-FFB
#'
#' @description
#' Infers a heuristic partition of observed variables into exogenous (\eqn{Y_2})
#' and endogenous (\eqn{Y_1}) blocks for use in NMF-FFB.
#' The method is based on positive-SEM logic, causal ordering, and optional
#' sign alignment using the first principal component (PC1).
#'
#' The procedure:
#' \itemize{
#'   \item internally standardizes variables (mean 0, sd 1),
#'   \item optionally flips signs so that most variables align positively with PC1,
#'   \item infers a causal ordering by repeatedly regressing each variable on the
#'         remaining ones and selecting the variable with the largest minimum
#'         standardized coefficient,
#'   \item determines an exogenous block by scanning the ordering from upstream
#'         and stopping at the first variable whose strongest parent coefficient
#'         exceeds \code{threshold}.
#' }
#'
#' If \code{n.exogenous} is supplied, it overrides the automatic threshold rule.
#'
#' @param x A numeric matrix or data frame with
#'   \strong{rows = samples} and \strong{columns = observed variables}.
#' @param n.exogenous Optional integer specifying the number of exogenous variables
#'   (\eqn{Y_2}). If \code{NULL}, the number is inferred automatically by the
#'   coefficient cut-off rule.
#' @param threshold Standardized regression-coefficient threshold used in the
#'   automatic exogenous–endogenous split. A variable is treated as endogenous
#'   once its maximum standardized parent coefficient exceeds this value.
#'   (Default: \code{0.1})
#' @param auto.flipped Logical; if \code{TRUE}, applies PC1-based automatic
#'   sign flipping after standardization to ensure consistent orientation.
#'   (Default: \code{TRUE})
#' @param verbose Logical; if \code{TRUE}, prints progress messages and the
#'   resulting variable split. (Default: \code{FALSE})
#' @param ... Reserved for future use; currently unused (also accepted
#'   by the \code{\link{nmf.ffb.split}} alias for argument forwarding).
#'
#' @return A list with:
#'   \item{endogenous.variables}{
#'     Character vector of variables selected as endogenous (\eqn{Y_1}).}
#'   \item{exogenous.variables}{
#'     Character vector of variables selected as exogenous (\eqn{Y_2}).}
#'   \item{ordered.variables}{
#'     Variables in inferred causal order (from exogenous to endogenous).}
#'   \item{is.flipped}{
#'     Logical vector indicating which variables were sign-flipped during processing.}
#'   \item{n.exogenous}{
#'     Integer giving the number of exogenous variables.}
#'
#' @examples
#' # Infer exogenous/endogenous split from iris
#' sp <- nmf.ffb.split(iris[, -5], n.exogenous = 2)
#' sp$endogenous.variables
#' sp$exogenous.variables
#'
#' @seealso \code{\link{nmf.ffb}}
#' @export
nmf.ffb.split <- function(x, n.exogenous = NULL, threshold = 0.1,
                          auto.flipped = TRUE, verbose = FALSE, ...) {

  if (!is.matrix(x) && !is.data.frame(x))
    stop("x must be a numeric matrix or data frame.")

  X_raw <- as.matrix(x)
  P <- ncol(X_raw)
  col_names <- colnames(X_raw)
  if (is.null(col_names)) {
    col_names <- paste0("V", 1:P)
    colnames(X_raw) <- col_names
  }

  # --------------------------------------------------------------------
  # Preprocessing Step 1: Standardize all variables
  #
  # Variables are centered and scaled (mean 0, sd 1). NMF-FFB requires
  # non-negative matrices, but the purpose of this function is only to
  # infer variable roles (Y1/Y2), so standardized values are allowed.
  #
  # Missing or NaN values resulting from constant columns are set to 0.
  # --------------------------------------------------------------------
  X_calc <- scale(X_raw, center = TRUE, scale = TRUE)
  X_calc[is.na(X_calc)] <- 0
  X_calc[is.nan(X_calc)] <- 0

  all_indices <- 1:P

  # --------------------------------------------------------------------
  # Preprocessing Step 2: Optional sign flipping based on PC1 alignment
  #
  # In positive SEM (and NMF-FFB), variables should ideally have
  # consistent sign orientation. To enforce this heuristic, variables
  # negatively correlated with the first principal component are flipped.
  #
  # This stabilizes the causal-ordering heuristic by avoiding mixtures
  # of arbitrary sign conventions in the raw data.
  # --------------------------------------------------------------------
  is.flipped <- rep(FALSE, P)
  names(is.flipped) <- col_names

  if (auto.flipped) {
    if (verbose) message("Step 0: Checking correlations with PC1 (on standardized data)...")

    svd_res <- svd(X_calc)
    pc1 <- svd_res$u[, 1]

    cors <- stats::cor(X_calc, pc1)
    # Ensure majority alignment with PC1
    if (stats::median(cors, na.rm = TRUE) < 0) cors <- -cors

    flip_idx <- which(cors < 0)

    if (length(flip_idx) > 0) {
      is.flipped[flip_idx] <- TRUE
      X_calc[, flip_idx] <- -X_calc[, flip_idx]

      if (verbose) {
        message(sprintf("   -> Detected %d flipped variables: %s",
                    length(flip_idx), paste(col_names[flip_idx], collapse=", ")))
      }
    }
  }

  # --------------------------------------------------------------------
  # Step 1: Causal ordering heuristic
  #
  # We infer an ordering of variables consistent with positive-SEM logic:
  # repeatedly select the variable that has the *largest minimum* coefficient
  # when regressed on the remaining variables. This favors variables that
  # are least explained by others → likely exogenous.
  #
  # The resulting order approximates a causal topological order in which
  # exogenous variables appear early and endogenous variables later.
  # --------------------------------------------------------------------
  if (verbose) message("Step 1: Inferring Causal Ordering...")

  active_set <- all_indices
  ordering_reversed <- integer(P)

  for (t in 1:(P - 1)) {
    scores <- numeric(length(active_set))

    for (i in seq_along(active_set)) {
      target_col <- active_set[i]
      pred_cols <- active_set[-i]

      y_vec <- X_calc[, target_col]
      X_mat <- X_calc[, pred_cols, drop = FALSE]

      coefs <- tryCatch({
        stats::coef(stats::lm(y_vec ~ X_mat - 1))
      }, error = function(e) rep(NA_real_, length(pred_cols)))

      # Positive SEM → keep the smallest coefficient (weakest positive predictor)
      if (all(is.na(coefs))) {
        scores[i] <- -Inf
      } else {
        scores[i] <- min(coefs, na.rm = TRUE)
      }
    }

    best_idx <- which.max(scores)
    ordering_reversed[t] <- active_set[best_idx]
    active_set <- active_set[-best_idx]
  }
  ordering_reversed[P] <- active_set[1]

  # Causal order: exogenous → endogenous
  ordering_indices <- rev(ordering_reversed)

  # --------------------------------------------------------------------
  # Step 2: Automatic identification of exogenous variables
  #
  # Sweep through the causal ordering. For each variable, regress it on all
  # earlier variables. If its strongest parent coefficient exceeds the
  # threshold, the variable is considered endogenous.
  #
  # Variables before this point → exogenous (Y2)
  # Variables after this point → endogenous (Y1)
  #
  # If n.exogenous is given, it overrides this automatic rule.
  # --------------------------------------------------------------------
  if (is.null(n.exogenous)) {
    if (verbose) message("Step 2: Detecting optimal cut-off for exogenous variables...")
    cutoff <- 1

    for (k in seq_len(max(P - 2, 0)) + 1) {
      curr_idx <- ordering_indices[k]
      parent_indices <- ordering_indices[seq_len(k - 1)]

      y_vec <- X_calc[, curr_idx]
      if (length(parent_indices) == 0) next
      X_parents <- X_calc[, parent_indices, drop = FALSE]

      coefs <- stats::coef(stats::lm(y_vec ~ X_parents - 1))
      coefs <- coefs[is.finite(coefs)]
      max_influence <- if (length(coefs) > 0) max(coefs) else 0

      if (max_influence > threshold) {
        if (verbose)
          message(sprintf("   [%d] %s : Max std.coef=%.3f -> Endogenous (Stop)",
                      k, col_names[curr_idx], max_influence))
        break
      } else {
        cutoff <- k
        if (verbose)
          message(sprintf("   [%d] %s : Max std.coef=%.3f -> Exogenous (Continue)",
                      k, col_names[curr_idx], max_influence))
      }
    }
    n.exogenous <- cutoff
  }

  # --------------------------------------------------------------------
  # Step 3: Final classification into Y1 and Y2
  #
  # Variables appearing early in the ordering (determined by cut-off) are
  # treated as exogenous (Y2). The remainder are endogenous (Y1).
  #
  # Ordered list shows the full inferred causal sequence.
  # --------------------------------------------------------------------
  if(n.exogenous < 1 || n.exogenous >= P)
    stop("n.exogenous must be between 1 and P-1.")
  idx_exo <- ordering_indices[1:n.exogenous]
  idx_endo <- ordering_indices[(n.exogenous + 1):P]

  exogenous.variables <- col_names[idx_exo]
  endogenous.variables <- col_names[idx_endo]
  ordered.variables <- col_names[ordering_indices]

  if (verbose) {
    message("\n--- Auto Split Result ---")
    message(sprintf("Exogenous (Y2, n=%d): %s",
                n.exogenous, paste(exogenous.variables, collapse=", ")))
    message(sprintf("Endogenous (Y1, n=%d): %s ...",
                length(endogenous.variables),
                paste(utils::head(endogenous.variables, 3), collapse=", ")))
  }

  return(list(
    endogenous.variables = endogenous.variables,
    exogenous.variables = exogenous.variables,
    ordered.variables = ordered.variables,
    is.flipped = is.flipped,
    n.exogenous = as.integer(n.exogenous)
  ))
}




############################################################
## Common DOT Helpers
############################################################

############################################################
## 1. nmf.ffb.DOT  (for NMF-FFB visualization)
############################################################

#' Generate a Graphviz DOT Diagram for an NMF-FFB Model
#'
#' @description
#' Creates a Graphviz DOT script that visualizes the structural network
#' estimated by \code{nmf.ffb}.
#' The resulting diagram displays:
#' \itemize{
#'   \item endogenous observed variables (\eqn{Y_1}),
#'   \item exogenous observed variables (\eqn{Y_2}),
#'   \item latent factors (\eqn{F_1}, \dots, \eqn{F_Q}),
#' }
#' together with the non-negative path coefficients whose magnitudes
#' exceed a user-specified threshold.
#'
#' Directed edges represent estimated relationships:
#' \itemize{
#'   \item \eqn{Y_2 \rightarrow F_q}: entries of \code{C2} (exogenous loadings),
#'   \item \eqn{F_q \rightarrow Y_1}: rows of \code{X} (factor-to-endogenous mappings),
#'   \item \eqn{Y_1 \rightarrow F_q}: entries of \code{C1} (feedback paths).
#' }
#'
#' Edge widths are scaled by coefficient magnitude, and nodes are placed
#' in optional visual clusters. Only variables participating in
#' edges above the threshold are displayed, while latent factors are always shown.
#'
#' @param result A list returned by \code{nmf.ffb}, containing matrices
#'   \code{X}, \code{C1}, and \code{C2}.
#' @param weight_scale Base scaling factor for edge widths.
#' @param weight_scale_c2 Scaling factor for edges
#'   \eqn{Y_2 \rightarrow F_q} (C2 matrix). Defaults to \code{weight_scale}.
#' @param weight_scale_x1 Scaling factor for edges
#'   \eqn{F_q \rightarrow Y_1} (X matrix). Defaults to \code{weight_scale}.
#' @param weight_scale_feedback Scaling factor for feedback edges
#'   \eqn{Y_1 \rightarrow F_q} (C1 matrix). Defaults to \code{weight_scale}.
#' @param threshold Minimum coefficient value needed for an edge to be drawn.
#' @param sig.level Significance level for filtering edges by p-value
#'   (requires inference results). Edges with p-value above this level are omitted.
#' @param rankdir Graphviz rank direction (e.g., \code{"LR"}, \code{"TB"}).
#' @param fill Logical; whether to use filled node shapes.
#' @param ... For backward compatibility: accepts deprecated names
#'   \code{weight_scale_y2f} (use \code{weight_scale_c2}) and
#'   \code{weight_scale_fy1} (use \code{weight_scale_x1}).
#' @param cluster.box Character string controlling the visibility and style
#'   of cluster frames around Y2, factors, and Y1 blocks.
#'   One of \code{"normal"}, \code{"faint"}, \code{"invisible"}, \code{"none"}.
#' @param cluster.labels Optional character vector of length 3 giving custom
#'   labels for the Y2, factor, and Y1 clusters.
#' @param hide.isolated Logical. If \code{TRUE} (default), Y1 and Y2 nodes
#'   that have no edges at or above \code{threshold} are excluded from the graph.
#' @param model Which fit of a \code{method = "fiml"} object to draw:
#'   \code{"selected"} (default; the BIC-selected model, with significance
#'   stars when inference results are present), \code{"null"} (the
#'   feed-forward null: \eqn{\Theta_1 = 0}, \eqn{\Theta_2} from
#'   \code{result$null}) or \code{"full"} (the unpenalized feedback fit from
#'   \code{result$full}).  The latter two ignore \code{result$coefficients},
#'   which refer to the selected model.  Ignored for \code{method = "mu"}
#'   objects, which carry a single fit.
#' @param sig.level Significance level for filtering structural edges
#'   (\eqn{C_1} feedback and \eqn{C_2} exogenous loadings) when
#'   inference results are present.  If \code{result} contains a
#'   \code{coefficients} data frame from \code{\link{nmf.ffb.inference}},
#'   only edges with \code{prob.unsupported < sig.level} are drawn, with
#'   significance stars (\code{*} \code{**} \code{***}) appended to
#'   the edge label.  The \eqn{X} (factor-to-\eqn{Y_1}) edges are
#'   never starred since the basis is not the inference target.
#'   Set to \code{NULL} to disable significance filtering and fall
#'   back to the \code{threshold} magnitude filter for both \eqn{C_1}
#'   and \eqn{C_2}.  Default is \code{0.1}.
#'
#' @return A character string representing a valid Graphviz DOT script.
#'
#' @examples
#' Y <- t(iris[, -5])
#' Y1 <- Y[1:2, ]
#' Y2 <- Y[3:4, ]
#' result <- nmf.ffb(Y1, Y2, rank = 2, maxit = 500)
#' dot <- nmf.ffb.DOT(result)
#' cat(dot)
#'
#' @seealso \code{\link{nmf.ffb}}, \code{\link{nmf.ffb.inference}},
#'   \code{\link{plot.nmfkc.DOT}}
#' @export
nmf.ffb.DOT <- function(result,
                        weight_scale          = 5,
                        weight_scale_c2       = weight_scale,
                        weight_scale_x1       = weight_scale,
                        weight_scale_feedback = weight_scale,
                        threshold             = 0.01,
                        sig.level             = 0.1,
                        rankdir               = "LR",
                        fill                  = TRUE,
                        cluster.box           = c("normal", "faint", "invisible", "none"),
                        cluster.labels        = NULL,
                        hide.isolated         = TRUE,
                        model                 = c("selected", "null", "full"),
                        ...) {

  ## Backward compatibility: accept deprecated names via ...
  extra_args <- base::list(...)
  if (!base::is.null(extra_args$weight_scale_y2f)) weight_scale_c2 <- extra_args$weight_scale_y2f
  if (!base::is.null(extra_args$weight_scale_fy1)) weight_scale_x1 <- extra_args$weight_scale_fy1

  ## Which fit of a fiml object to draw.  The null and the unpenalized fit
  ## are substituted into C1 / C2; the inference table (selected model) is
  ## dropped so that no stars from another model decorate their edges.
  model <- match.arg(model)
  if (model != "selected") {
    alt <- result[[model]]
    if (is.null(alt))
      stop(sprintf("result has no `%s` component; model = \"%s\" needs an nmf.ffb(method = \"fiml\") fit.",
                   model, model))
    result$C1 <- if (model == "null") 0 * result$C1 else alt$C1
    result$C2 <- alt$C2
    result$coefficients <- NULL
  }

  ## ---------------------------------------------------------------
  ## Cluster style selection
  ## ---------------------------------------------------------------
  cluster.box <- match.arg(cluster.box)

  cluster_style <- switch(cluster.box,
                          normal    = "rounded",
                          faint     = "rounded,dashed",
                          invisible = "rounded",
                          none      = "none")

  cluster_color <- switch(cluster.box,
                          normal    = "black",
                          faint     = "gray80",
                          invisible = "none",
                          none      = "none")

  cluster_penwidth <- switch(cluster.box,
                             normal    = 1.0,
                             faint     = 0.7,
                             invisible = 0.0,
                             none      = 0.0)

  ## ---------------------------------------------------------------
  ## Cluster titles
  ## ---------------------------------------------------------------
  default.labels <- c("Exogenous (Y2)", "Latent Factors", "Endogenous (Y1)")
  if (is.null(cluster.labels)) {
    titles <- default.labels
  } else {
    if (!is.character(cluster.labels) || length(cluster.labels) != 3L) {
      stop("cluster.labels must be a character vector of length 3: c(label_Y2, label_F, label_Y1).")
    }
    titles <- cluster.labels
  }

  ## ---------------------------------------------------------------
  ## Extract matrices and labels
  ## ---------------------------------------------------------------
  X  <- result$X
  C1 <- result$C1
  C2 <- result$C2

  if (is.null(X) || is.null(C1) || is.null(C2)) {
    stop("result must contain elements X, C1, and C2.")
  }

  Y1_labels <- rownames(X)
  Y2_labels <- colnames(C2)

  if (is.null(Y1_labels)) Y1_labels <- paste0("Y1_", seq_len(nrow(X)))
  if (is.null(Y2_labels)) Y2_labels <- paste0("Y2_", seq_len(ncol(C2)))

  P1 <- length(Y1_labels)
  P2 <- length(Y2_labels)
  Q  <- ncol(X)

  ## ---------------------------------------------------------------
  ## Identify nodes involved in edges >= threshold
  ## ---------------------------------------------------------------
  if (isTRUE(hide.isolated)) {
    used_y2 <- apply(C2, 2L, function(col) any(col >= threshold, na.rm = TRUE))
    used_y1_from_X  <- apply(X,  1L, function(row) any(row >= threshold, na.rm = TRUE))
    used_y1_from_C1 <- apply(C1, 2L, function(col) any(col >= threshold, na.rm = TRUE))
    used_y1 <- used_y1_from_X | used_y1_from_C1
    idx_y1 <- which(used_y1)
    idx_y2 <- which(used_y2)
  } else {
    idx_y1 <- seq_len(P1)
    idx_y2 <- seq_len(P2)
  }

  ## ---------------------------------------------------------------
  ## Assign internal DOT node IDs
  ## ---------------------------------------------------------------
  Y1_ids <- .nmfkc_dot_sanitize_id(paste0("Y1_", seq_len(P1)))
  Y2_ids <- .nmfkc_dot_sanitize_id(paste0("Y2_", seq_len(P2)))
  F_ids  <- .nmfkc_dot_sanitize_id(paste0("F_",  seq_len(Q)))

  ## Node colors
  COLOR_Y2_NODE     <- "lightcoral"
  COLOR_Y1_NODE     <- "lightblue"
  COLOR_FACTOR_NODE <- "wheat"

  ## ---------------------------------------------------------------
  ## Header
  ## ---------------------------------------------------------------
  dot_script <- .nmfkc_dot_header(
    graph_name = "NMF_SEM_Full_Mechanism",
    rankdir    = rankdir
  )

  ## ---------------------------------------------------------------
  ## Y2 cluster
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    "\n  // Exogenous variables (Y2)\n",
    .nmfkc_dot_cluster_nodes(
      cluster_id        = "Y2",
      title             = titles[1],
      node_ids          = Y2_ids[idx_y2],
      node_labels       = Y2_labels[idx_y2],
      shape             = "box",
      fill              = fill,
      fillcolor         = COLOR_Y2_NODE,
      line_width        = 1.5,
      cluster_style     = cluster_style,
      cluster_color     = cluster_color,
      cluster_penwidth  = cluster_penwidth
    )
  )

  ## ---------------------------------------------------------------
  ## Y1 cluster
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    "\n  // Endogenous variables (Y1)\n",
    .nmfkc_dot_cluster_nodes(
      cluster_id        = "Y1",
      title             = titles[3],
      node_ids          = Y1_ids[idx_y1],
      node_labels       = Y1_labels[idx_y1],
      shape             = "box",
      fill              = fill,
      fillcolor         = COLOR_Y1_NODE,
      line_width        = 1.5,
      cluster_style     = cluster_style,
      cluster_color     = cluster_color,
      cluster_penwidth  = cluster_penwidth
    )
  )

  ## ---------------------------------------------------------------
  ## Latent factor cluster
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    "\n  // Latent Factors (F)\n",
    .nmfkc_dot_cluster_nodes(
      cluster_id        = "F",
      title             = titles[2],
      node_ids          = F_ids,
      node_labels       = paste0("Factor ", seq_len(Q)),
      shape             = "ellipse",
      fill              = fill,
      fillcolor         = COLOR_FACTOR_NODE,
      line_width        = 1.0,
      cluster_style     = cluster_style,
      cluster_color     = cluster_color,
      cluster_penwidth  = cluster_penwidth
    )
  )

  ## ---------------------------------------------------------------
  ## Edge defaults
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    '\n  edge [fontname="Arial", fontsize=8, arrowhead=open];\n'
  )

  pw     <- .nmfkc_dot_penwidth
  digits <- .nmfkc_dot_digits_from_threshold(threshold)
  fmtc   <- function(x) .nmfkc_dot_format_coef(x, digits)

  ## ---------------------------------------------------------------
  ## Significance stars for C1 (feedback) and C2 (exogenous) edges,
  ## both from nmf.ffb.inference().  X (F -> Y1) edges are NOT
  ## starred even when inference results are present, since the
  ## basis is not the inference target.
  ##
  ## Filters by the optional `Type` column ("C1" / "C2") if present
  ## (newer inference output); falls back to row-name matching for
  ## back-compatibility with pre-v0.6.8 inference output that only
  ## carried C2.
  ## ---------------------------------------------------------------
  pval_to_stars <- function(p) {
    if (!is.finite(p)) ""
    else if (p < 0.001) "***"
    else if (p < 0.01)  "**"
    else if (p < 0.05)  "*"
    else                ""
  }

  C1_stars <- NULL; C1_show <- NULL
  C2_stars <- NULL; C2_show <- NULL
  if (!is.null(result$coefficients)) {
    cf        <- result$coefficients
    fac_names <- rownames(C2)
    y1_names  <- colnames(C1)
    y2_names  <- colnames(C2)

    has_type <- !is.null(cf$Type)
    cf_C1 <- if (has_type) cf[cf$Type == "C1", , drop = FALSE] else NULL
    cf_C2 <- if (has_type) cf[cf$Type == "C2", , drop = FALSE] else cf

    ## ---- C2 stars (Y2 -> F) ----
    C2_stars <- matrix("", nrow = Q, ncol = P2)
    C2_pval  <- matrix(NA_real_, nrow = Q, ncol = P2)
    for (k in seq_len(nrow(cf_C2))) {
      q  <- match(cf_C2$Basis[k], fac_names)
      p2 <- match(cf_C2$Covariate[k], y2_names)
      if (!is.na(q) && !is.na(p2) && !is.na(cf_C2$prob.unsupported[k])) {
        C2_pval[q, p2]  <- cf_C2$prob.unsupported[k]
        C2_stars[q, p2] <- pval_to_stars(cf_C2$prob.unsupported[k])
      }
    }
    if (!is.null(sig.level)) {
      C2_show <- !is.na(C2_pval) & C2_pval < sig.level
    }

    ## ---- C1 stars (Y1 -> F, feedback) ----
    if (!is.null(cf_C1) && nrow(cf_C1) > 0L) {
      C1_stars <- matrix("", nrow = Q, ncol = P1)
      C1_pval  <- matrix(NA_real_, nrow = Q, ncol = P1)
      for (k in seq_len(nrow(cf_C1))) {
        q  <- match(cf_C1$Basis[k], fac_names)
        p1 <- match(cf_C1$Covariate[k], y1_names)
        if (!is.na(q) && !is.na(p1) && !is.na(cf_C1$prob.unsupported[k])) {
          C1_pval[q, p1]  <- cf_C1$prob.unsupported[k]
          C1_stars[q, p1] <- pval_to_stars(cf_C1$prob.unsupported[k])
        }
      }
      if (!is.null(sig.level)) {
        C1_show <- !is.na(C1_pval) & C1_pval < sig.level
      }
    }
  }

  ## ---------------------------------------------------------------
  ## 1. Y2 → F edges (C2)
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    '\n  // 1. External Driving (Y2 -> Factor) [C2]\n',
    '  edge [color=black, fontcolor=black, style=solid];\n'
  )

  max_C2 <- suppressWarnings(max(C2, na.rm = TRUE))
  if (is.finite(max_C2) && max_C2 > 0) {
    for (q in seq_len(Q)) {
      for (p2 in idx_y2) {
        weight <- C2[q, p2]
        show <- if (!is.null(C2_show)) C2_show[q, p2]
                else is.finite(weight) && weight >= threshold
        if (show) {
          pen <- pw(weight, max_C2, weight_scale_c2)
          lab <- fmtc(weight)
          if (!is.null(C2_stars)) lab <- paste0(lab, C2_stars[q, p2])
          path <- sprintf('  %s -> %s [label="%s", penwidth=%.2f];\n',
                          Y2_ids[p2], F_ids[q], lab, pen)
          dot_script <- paste0(dot_script, path)
        }
      }
    }
  }

  ## ---------------------------------------------------------------
  ## 2. F → Y1 edges (X)
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    '\n  // 2. Generation (Factor -> Y1) [X]\n',
    '  edge [color="gray0", fontcolor="gray0", style=solid];\n'
  )

  max_X <- suppressWarnings(max(X, na.rm = TRUE))
  if (is.finite(max_X) && max_X > 0) {
    for (q in seq_len(Q)) {
      for (p1 in idx_y1) {
        weight <- X[p1, q]
        if (is.finite(weight) && weight >= threshold) {
          pen <- pw(weight, max_X, weight_scale_x1)
          lab <- fmtc(weight)
          path <- sprintf('  %s -> %s [label="%s", penwidth=%.2f];\n',
                          F_ids[q], Y1_ids[p1], lab, pen)
          dot_script <- paste0(dot_script, path)
        }
      }
    }
  }

  ## ---------------------------------------------------------------
  ## 3. Y1 → F edges (C1 feedback)
  ## ---------------------------------------------------------------
  dot_script <- paste0(
    dot_script,
    '\n  // 3. Internal Feedback (Y1 -> Factor) [C1]\n',
    '  edge [style=dashed, color="gray0", fontcolor="gray0"];\n'
  )

  max_C1 <- suppressWarnings(max(C1, na.rm = TRUE))
  if (is.finite(max_C1) && max_C1 > 0) {
    for (q in seq_len(Q)) {
      for (p1 in idx_y1) {
        weight <- C1[q, p1]
        ## When inference results provide C1 p-values, prefer the
        ## significance-based filter (parallels the C2 branch above);
        ## otherwise fall back to the magnitude threshold.
        show <- if (!is.null(C1_show)) C1_show[q, p1]
                else is.finite(weight) && weight >= threshold
        if (show) {
          pen <- pw(weight, max_C1, weight_scale_feedback)
          lab <- fmtc(weight)
          if (!is.null(C1_stars)) lab <- paste0(lab, C1_stars[q, p1])
          path <- sprintf('  %s -> %s [label="%s", penwidth=%.2f];\n',
                          Y1_ids[p1], F_ids[q], lab, pen)
          dot_script <- paste0(dot_script, path)
        }
      }
    }
  }

  result <- paste0(dot_script, "}\n")
  class(result) <- c("nmf.ffb.DOT", "nmfkc.DOT")
  result
}





############################################################
## 2. nmfkc.DOT  (Static NMF / NMF-with-covariates visualization)
############################################################

#' Generate Graphviz DOT Scripts for NMF or NMF-with-Covariates Models
#'
#' @description
#' Produces a Graphviz DOT script visualizing the structure of an NMF model
#' (\eqn{Y \approx X C A}) or its simplified forms.
#'
#' Supported visualization types:
#' \itemize{
#'   \item \code{"YX"} — Standard NMF view: latent factors \eqn{X} map to observations \eqn{Y}.
#'   \item \code{"YA"} — Direct regression view: covariates \eqn{A} map directly to \eqn{Y}
#'         using the combined coefficient matrix \eqn{X C}.
#'   \item \code{"YXA"} — Full tri-factorization: \eqn{A \rightarrow C \rightarrow X \rightarrow Y}.
#' }
#'
#' Edge widths are scaled by coefficient magnitude, and nodes with no edges
#' above the threshold are omitted from the visualization.
#'
#' @param result The return value from \code{nmfkc}, containing matrices
#'   \code{X}, \code{B}, and optionally \code{C}.
#' @param type Character string specifying the visualization style:
#'   one of \code{"YX"}, \code{"YA"}, \code{"YXA"}.
#' @param threshold Minimum coefficient magnitude to display an edge. When
#'   \code{C.signed = TRUE} (signed \eqn{\Theta}) this is an absolute-value cut,
#'   i.e. an edge is drawn when \eqn{|coef| \ge} \code{threshold}.
#' @param C.signed Logical or \code{NULL}. Whether the coefficient matrix
#'   \eqn{\Theta} (\code{= C}) may be signed (real-valued). When \code{TRUE},
#'   the \code{threshold} is applied to \eqn{|coef|}, edge widths are scaled by
#'   \eqn{|coef|}, and negative edges are drawn as black \strong{dashed} lines
#'   (positive edges remain solid) to distinguish them; the numeric edge labels
#'   keep their sign. When \code{FALSE}, the historical non-negative behaviour is used
#'   (negative coefficients fall below the threshold and are hidden). When
#'   \code{NULL} (default), the mode is auto-detected from the fit
#'   (\code{result$C.signed}) or from the presence of negative entries in
#'   \eqn{C} / \eqn{XC}. The basis \eqn{X} is always non-negative, so
#'   \eqn{X \rightarrow Y} edges are unaffected.
#' @param sig.level Significance level for filtering C edges when inference
#'   results are available (i.e., \code{x} is of class \code{"nmfkc.inference"}).
#'   Only edges with p-value below \code{sig.level} are shown, decorated with
#'   significance stars (\code{*}, \code{**}, \code{***}). Set to \code{NULL}
#'   to disable filtering and show all edges above \code{threshold}. Default is 0.1.
#' @param rankdir Graphviz rank direction (e.g., \code{"LR"}, \code{"TB"}).
#' @param fill Logical; whether nodes should be drawn with filled shapes.
#' @param weight_scale Base scaling factor for edge widths.
#' @param weight_scale_ax Scaling factor for edges \eqn{A \rightarrow X} (type \code{"YXA"}).
#' @param weight_scale_xy Scaling factor for edges \eqn{X \rightarrow Y}.
#' @param weight_scale_ay Scaling factor for edges \eqn{A \rightarrow Y} (type \code{"YA"}).
#' @param Y.label Optional character vector for labels of Y nodes.
#' @param X.label Optional character vector for labels of X (latent factor) nodes.
#' @param A.label Optional character vector for labels of A (covariate) nodes.
#' @param Y.title Cluster title for Y nodes.
#' @param X.title Cluster title for X nodes.
#' @param A.title Cluster title for A nodes.
#' @param hide.isolated Logical. If \code{TRUE} (default), Y and A nodes that have no
#'   edges at or above \code{threshold} are excluded from the graph.
#'
#' @return A character string representing a Graphviz DOT script.
#'
#' @seealso \code{\link{nmfkc}}, \code{\link{nmf.rrr.DOT}}, \code{\link{nmf.ffb.DOT}},
#'   \code{\link{nmfkc.ar.DOT}}, \code{\link{plot.nmfkc.DOT}}
#' @examples
#' Y <- matrix(cars$dist, nrow = 1)
#' A <- rbind(1, cars$speed)
#' result <- nmfkc(Y, A, rank = 1)
#' dot <- nmfkc.DOT(result)
#' cat(dot)
#'
#' @export
nmfkc.DOT <- function(
    result,
    type = c("YX","YA","YXA"),
    threshold = 0.01,
    C.signed = NULL,
    sig.level = 0.1,
    rankdir   = "LR",
    fill      = TRUE,
    weight_scale    = 5,
    weight_scale_ax = weight_scale,
    weight_scale_xy = weight_scale,
    weight_scale_ay = weight_scale,
    Y.label = NULL, X.label = NULL, A.label = NULL,
    Y.title = "Observation (Y)",
    X.title = "Basis (X)",
    A.title = "Covariates (A)",
    hide.isolated = TRUE
) {

  type <- match.arg(type)

  ## ---------------------------------------------------------
  ## Required matrices
  ## ---------------------------------------------------------
  X <- result$X
  B <- result$B
  if (is.null(X) || is.null(B)) {
    stop("result must contain X and B.")
  }

  ## If C exists and is a proper NMF-with-covariates factor:
  hasA <- !is.null(result$C) && ncol(result$C) != ncol(B)

  P <- nrow(X)
  Q <- ncol(X)

  ## ---------------------------------------------------------
  ## Labels
  ## ---------------------------------------------------------
  Y_labels <- if (is.null(Y.label)) rownames(X) else Y.label
  X_labels <- if (is.null(X.label)) colnames(X) else X.label

  if (is.null(Y_labels)) Y_labels <- paste0("Y", seq_len(P))
  if (is.null(X_labels)) X_labels <- paste0("Factor", seq_len(Q))

  Y_ids <- .nmfkc_dot_sanitize_id(paste0("Y_", seq_len(P)))
  X_ids <- .nmfkc_dot_sanitize_id(paste0("X_", seq_len(Q)))

  ## ---------------------------------------------------------
  ## Covariates and tri-factorization
  ## ---------------------------------------------------------
  if (hasA) {
    C <- as.matrix(result$C)          # Q x R
    A_cols <- ncol(C)
    A_labels <- if (is.null(A.label)) colnames(C) else A.label
    if (is.null(A_labels)) A_labels <- paste0("A", seq_len(A_cols))
    A_ids <- .nmfkc_dot_sanitize_id(paste0("A_", seq_len(A_cols)))

    ## Combined mapping for type = "YA"
    XC_mat <- X %*% C   # P x R

  } else if (type == "YX") {
    ## No A block needed
    C <- NULL
    A_cols <- 0L
    A_labels <- NULL
    A_ids <- NULL
    XC_mat <- NULL

  } else {
    stop("The model structure (matrix C) is incompatible with type 'YXA' or 'YA'.")
  }

  ## ---------------------------------------------------------
  ## Sign mode of C (= Theta). When NULL, auto-detect from the fit
  ## (result$C.signed) or from negative entries in C / XC. In signed mode the
  ## threshold is an absolute-value cut (|coef| >= threshold), edge widths use
  ## |coef|, and negative edges are coloured to distinguish them. The basis X
  ## is always non-negative, so X -> Y edges are unaffected.
  ## ---------------------------------------------------------
  if (is.null(C.signed)) {
    C.signed <- isTRUE(result$C.signed) ||
      (!is.null(C)      && any(C      < 0, na.rm = TRUE)) ||
      (!is.null(XC_mat) && any(XC_mat < 0, na.rm = TRUE))
  }
  C.signed <- isTRUE(C.signed)
  ## magnitude used for thresholding and penwidth scaling (|.| in signed mode)
  mag <- if (C.signed) function(v) abs(v) else function(v) v
  ## edge line style by sign (negative -> dashed) in signed mode; colour stays black
  estyle <- function(v) if (C.signed && is.finite(v) && v < 0) "dashed" else "solid"

  ## ---------------------------------------------------------
  ## Filter isolated nodes (hide.isolated)
  ## ---------------------------------------------------------
  idx_Y <- seq_len(P)
  idx_A <- if (hasA) seq_len(A_cols) else integer(0)
  if (isTRUE(hide.isolated)) {
    if (type == "YX" || type == "YXA") {
      # Y is connected via X: X[i, ] >= threshold
      used_Y <- apply(X, 1L, function(row) any(row >= threshold, na.rm = TRUE))
      idx_Y <- which(used_Y)
    } else if (type == "YA" && !is.null(XC_mat)) {
      # Y is connected via XC: |XC_mat[i, ]| >= threshold
      used_Y <- apply(XC_mat, 1L, function(row) any(mag(row) >= threshold, na.rm = TRUE))
      idx_Y <- which(used_Y)
    }
    if (hasA && type %in% c("YXA", "YA")) {
      if (type == "YXA" && !is.null(C)) {
        # A is connected via C: |C[, k]| >= threshold
        used_A <- apply(C, 2L, function(col) any(mag(col) >= threshold, na.rm = TRUE))
        idx_A <- which(used_A)
      } else if (type == "YA" && !is.null(XC_mat)) {
        # A is connected via XC: |XC_mat[, k]| >= threshold
        used_A <- apply(XC_mat, 2L, function(col) any(mag(col) >= threshold, na.rm = TRUE))
        idx_A <- which(used_A)
      }
    }
  }

  ## ---------------------------------------------------------
  ## DOT header
  ## ---------------------------------------------------------
  scr <- .nmfkc_dot_header(graph_name = "NMF", rankdir = rankdir)

  ## ---------------------------------------------------------
  ## Y node cluster
  ## ---------------------------------------------------------
  scr <- paste0(
    scr,
    '\n  // Output variables (Y)\n',
    .nmfkc_dot_cluster_nodes(
      cluster_id  = "Y",
      title       = Y.title,
      node_ids    = Y_ids[idx_Y],
      node_labels = Y_labels[idx_Y],
      shape       = "box",
      fill        = fill,
      fillcolor   = "lightblue",
      line_width  = 1.5
    )
  )

  ## ---------------------------------------------------------
  ## X node cluster
  ## (Hidden when type = "YA")
  ## ---------------------------------------------------------
  if (type != "YA") {
    scr <- paste0(
      scr,
      '\n  // Latent factors (X)\n',
      .nmfkc_dot_cluster_nodes(
        cluster_id  = "X",
        title       = X.title,
        node_ids    = X_ids,
        node_labels = X_labels,
        shape       = "ellipse",
        fill        = fill,
        fillcolor   = "wheat",
        line_width  = 1.0
      )
    )
  }

  ## ---------------------------------------------------------
  ## A node cluster
  ## (Only for "YA" and "YXA")
  ## ---------------------------------------------------------
  if (type != "YX" && hasA) {
    scr <- paste0(
      scr,
      '\n  // Covariates (A)\n',
      .nmfkc_dot_cluster_nodes(
        cluster_id  = "A",
        title       = A.title,
        node_ids    = A_ids[idx_A],
        node_labels = A_labels[idx_A],
        shape       = "box",
        fill        = fill,
        fillcolor   = "lightcoral",
        line_width  = 1.5
      )
    )
  }

  ## ---------------------------------------------------------
  ## Significance stars for C edges (from nmfkc.inference)
  ## ---------------------------------------------------------
  C_stars <- NULL
  C_show  <- NULL
  if (!is.null(result$coefficients) && hasA) {
    C_stars <- matrix("", nrow = Q, ncol = A_cols)
    C_pval  <- matrix(NA_real_, nrow = Q, ncol = A_cols)
    cf <- result$coefficients
    basis_names <- rownames(C)
    cov_names   <- colnames(C)
    ## nmfkc.net.inference names the two factors Basis.row / Basis.col, because
    ## there both axes are bases.  Its result still inherits class "nmfkc", so it
    ## can reach here; without this fallback match() returned integer(0) and the
    ## `if` below errored with "argument is of length zero".  Same fallback that
    ## nmfkc.net.DOT already carries.
    row_col <- if ("Basis.row" %in% names(cf)) "Basis.row" else "Basis"
    col_col <- if ("Basis.col" %in% names(cf)) "Basis.col" else "Covariate"
    for (k in seq_len(nrow(cf))) {
      q <- match(cf[[row_col]][k], basis_names)
      r <- match(cf[[col_col]][k], cov_names)
      if (!is.na(q) && !is.na(r) && !is.na(cf$p_value[k])) {
        p <- cf$p_value[k]
        C_pval[q, r] <- p
        if (p < 0.001)      C_stars[q, r] <- "***"
        else if (p < 0.01)  C_stars[q, r] <- "**"
        else if (p < 0.05)  C_stars[q, r] <- "*"
      }
    }
    if (!is.null(sig.level)) {
      C_show <- !is.na(C_pval) & C_pval < sig.level
    }
  }

  ## ---------------------------------------------------------
  ## Edge defaults
  ## ---------------------------------------------------------
  scr <- paste0(
    scr,
    '\n  edge [fontname="Arial", fontsize=8, arrowhead=open];\n'
  )

  pw     <- .nmfkc_dot_penwidth
  digits <- .nmfkc_dot_digits_from_threshold(threshold)
  fmtc   <- function(x) .nmfkc_dot_format_coef(x, digits)


  ## =========================================================
  ## Case 1: Full tri-factorization (A → X → Y)
  ## =========================================================
  if (type == "YXA") {

    ## ---- A → X (C) edges ----
    scr <- paste0(
      scr,
      '\n  // A -> X edges (C)\n',
      '  edge [color=black, fontcolor=black, style=solid];\n'
    )

    max_C <- suppressWarnings(max(mag(C), na.rm = TRUE))
    if (is.finite(max_C) && max_C > 0) {
      for (q in seq_len(Q)) {
        for (k in idx_A) {
          val <- C[q, k]
          show <- if (!is.null(C_show)) C_show[q, k]
                  else is.finite(val) && mag(val) >= threshold
          if (show) {
            pen <- pw(mag(val), max_C, weight_scale_ax)
            lab <- fmtc(val)
            if (!is.null(C_stars)) lab <- paste0(lab, C_stars[q, k])
            sty <- estyle(val)
            scr <- paste0(
              scr,
              sprintf('  %s -> %s [label="%s", penwidth=%.2f, style=%s];\n',
                      A_ids[k], X_ids[q], lab, pen, sty)
            )
          }
        }
      }
    }

    ## ---- X → Y edges ----
    scr <- paste0(
      scr,
      '\n  // X -> Y edges (X)\n',
      '  edge [color="gray0", fontcolor="gray0", style=solid];\n'
    )

    max_X <- suppressWarnings(max(X, na.rm = TRUE))
    if (is.finite(max_X) && max_X > 0) {
      for (i in idx_Y) {
        for (j in seq_len(Q)) {
          val <- X[i, j]
          if (is.finite(val) && val >= threshold) {
            pen <- pw(val, max_X, weight_scale_xy)
            lab <- fmtc(val)
            scr <- paste0(
              scr,
              sprintf('  %s -> %s [label="%s", penwidth=%.2f];\n',
                      X_ids[j], Y_ids[i], lab, pen)
            )
          }
        }
      }
    }


    ## =========================================================
    ## Case 2: Direct regression view (A → Y)
    ## =========================================================
  } else if (type == "YA") {

    scr <- paste0(
      scr,
      '\n  // A -> Y edges (X %*% C)\n',
      '  edge [color=black, fontcolor=black, style=solid];\n'
    )

    max_XC <- suppressWarnings(max(mag(XC_mat), na.rm = TRUE))
    if (is.finite(max_XC) && max_XC > 0) {
      for (i in idx_Y) {
        for (k in idx_A) {
          val <- XC_mat[i, k]
          if (is.finite(val) && mag(val) >= threshold) {
            pen <- pw(mag(val), max_XC, weight_scale_ay)
            lab <- fmtc(val)
            sty <- estyle(val)
            scr <- paste0(
              scr,
              sprintf('  %s -> %s [label="%s", penwidth=%.2f, style=%s];\n',
                      A_ids[k], Y_ids[i], lab, pen, sty)
            )
          }
        }
      }
    }


    ## =========================================================
    ## Case 3: Standard NMF (X → Y)
    ## =========================================================
  } else if (type == "YX") {

    scr <- paste0(
      scr,
      '\n  // X -> Y edges (X)\n',
      '  edge [color="gray0", fontcolor="gray0", style=solid];\n'
    )

    max_X <- suppressWarnings(max(X, na.rm = TRUE))
    if (is.finite(max_X) && max_X > 0) {
      for (i in idx_Y) {
        for (j in seq_len(Q)) {
          val <- X[i, j]
          if (is.finite(val) && val >= threshold) {
            pen <- pw(val, max_X, weight_scale_xy)
            lab <- fmtc(val)
            scr <- paste0(
              scr,
              sprintf('  %s -> %s [label="%s", penwidth=%.2f];\n',
                      X_ids[j], Y_ids[i], lab, pen)
            )
          }
        }
      }
    }
  }

  result <- paste0(scr, "}\n")
  class(result) <- "nmfkc.DOT"
  result
}


#' @title Plot method for nmfkc.DOT objects
#' @description
#' Renders a DOT graph string using \code{DiagrammeR::grViz}.
#' If the \pkg{DiagrammeR} package is not installed, prints the DOT source
#' to the console instead.
#'
#' This method handles all DOT objects produced by the nmfkc package:
#' \code{\link{nmfkc.DOT}}, \code{\link{nmf.rrr.DOT}}, \code{\link{nmf.ffb.DOT}},
#' and \code{\link{nmfkc.ar.DOT}}.
#'
#' @param x An object of class \code{"nmfkc.DOT"} (or a subclass thereof).
#' @param ... Not used.
#' @return Called for its side effect (rendering). Returns \code{x} invisibly.
#' @seealso \code{\link{nmfkc.DOT}}, \code{\link{nmf.rrr.DOT}},
#'   \code{\link{nmf.ffb.DOT}}, \code{\link{nmfkc.ar.DOT}},
#'   \code{\link{print.nmfkc.DOT}}
#' @export
plot.nmfkc.DOT <- function(x, ...) {
  if (requireNamespace("DiagrammeR", quietly = TRUE)) {
    print(DiagrammeR::grViz(as.character(x)))
  } else {
    message("DiagrammeR package not installed. Printing DOT source:")
    cat(as.character(x))
  }
  invisible(x)
}


#' @title Print method for nmfkc.DOT objects
#' @description
#' Writes the Graphviz source out as text.  Without it, printing a DOT object
#' fell to \code{print.default}, which showed the whole graph as one escaped
#' string (\code{[1] "digraph NMF \{\\n ..."}) --- unreadable.  Shared by every
#' \code{*.DOT} function through the common \code{"nmfkc.DOT"} class.
#' @param x An object of class \code{"nmfkc.DOT"} (or a subclass).
#' @param ... Not used.
#' @return \code{x}, invisibly.
#' @seealso \code{\link{plot.nmfkc.DOT}}, \code{\link{nmfkc.DOT}}
#' @export
print.nmfkc.DOT <- function(x, ...) {
  cat(as.character(x))
  cat("# Use plot() to render (requires DiagrammeR),",
      "or write this text to a .dot file.\n")
  invisible(x)
}


#' @title Attach the identification diagnostic to an nmf.ffb.inference result (Internal)
#' @description \eqn{Y_1 \approx X(\Theta_1 Y_1 + \Theta_2 Y_2)} is
#'   \eqn{Y_1 \approx X\Theta A} with \eqn{\Theta = [\Theta_1\ \Theta_2]} and
#'   \eqn{A = [Y_1; Y_2]}, all of \eqn{X}, \eqn{\Theta_1}, \eqn{\Theta_2} non-negative
#'   and the columns of \eqn{X} summing to one, so the diagnostic of
#'   \code{nmfkc.inference()} applies: the bootstrap is conditional on \eqn{\hat X},
#'   and \code{identified.approx} says whether \eqn{(X, \Theta)} is unique up to a
#'   permutation.  The exclusion restriction's zeros in \eqn{\Theta_1} help.  The
#'   full covariance \eqn{\Phi} of the latent errors absorbs any \eqn{T}, so it
#'   adds nothing.  An algebraic check on the point estimate: no re-fits.
#' @param object The inference result.
#' @param Y1,Y2 The data.
#' @param ex The \code{...} list (reads \code{sep.tol}, \code{pin.tol}).
#' @keywords internal
#' @noRd
.ffb.attach.ident <- function(object, Y1, Y2, ex) {
  v <- .ident.verdict(object$X, base::cbind(object$C1, object$C2),
                      base::rbind(base::as.matrix(Y1), base::as.matrix(Y2)),
                      scale.fixed = TRUE,
                      sep.tol = if (!is.null(ex$sep.tol)) ex$sep.tol else 0.999,
                      pin.tol = if (!is.null(ex$pin.tol)) ex$pin.tol else 0.01)
  object$identified.approx <- v$identified.approx
  object$identification    <- v$identification
  object
}