/*──────────────────────────────────────────────────────────────────────────────
  trop_data_transfer.mata

  Data transfer layer between Mata and the native-code plugin.

  Stores variable indices, regularization grids, estimation options, and
  pre-allocated output matrices into the Stata namespace (scalars, matrices,
  global macros) so that the plugin can retrieve them via the Stata Plugin
  Interface (SPI).  Observation-level data is read by the plugin directly
  through SF_vdata(); this module does not construct large Mata matrices.

  Contents
    trop_prepare_data()              variable-index and panel-dimension setup
    trop_prepare_lambda_grids()      regularization-parameter grids
    trop_prepare_options()           estimation and inference options
    trop_prepare_output_matrices()   pre-allocate plugin output matrices
    trop_prepare_bootstrap()         bootstrap parameters and output matrix
    trop_cleanup_temp_vars()         drop all __trop_* temporaries
──────────────────────────────────────────────────────────────────────────────*/

version 17
mata:
mata set matastrict on

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_data()

  Maps variable names to positional indices within the plugin-call varlist
  and stores them as Stata scalars.  Records panel dimensions (N, T) and
  constructs the T x T time-distance matrix used by the exponential
  time-decay kernel  theta_s(t) = exp(-lambda_time * |t - s|).

  Arguments
    y_varname       outcome variable name          (Y)
    d_varname       binary treatment indicator name (W)
    panel_varname   panel (unit) identifier name
    time_varname    time period identifier name
    n_units         number of cross-sectional units (N)
    n_periods       number of time periods          (T)

  Stored scalars
    __trop_y_varindex, __trop_d_varindex, __trop_ctrl_varindex,
    __trop_panel_varindex, __trop_time_varindex,
    __trop_n_units, __trop_n_periods

  Stored globals
    __trop_varlist, __trop_y_varname, __trop_d_varname,
    __trop_panel_varname, __trop_time_varname

  Stored matrix
    __trop_time_dist   T x T absolute period distances
──────────────────────────────────────────────────────────────────────────────*/

void trop_prepare_data(
    string scalar y_varname,
    string scalar d_varname,
    string scalar panel_varname,
    string scalar time_varname,
    real scalar n_units,
    real scalar n_periods
)
{
    /* Varlist order for `plugin call`: Y, W, panel_id, time_id.
       Variable indices are therefore fixed at 1..4. */
    st_global("__trop_varlist", y_varname + " " + d_varname + " " + panel_varname + " " + time_varname)

    st_numscalar("__trop_y_varindex", 1)
    st_numscalar("__trop_d_varindex", 2)

    /* Control indicator reuses the treatment column; the plugin
       identifies control observations as those with W_{it} = 0. */
    st_numscalar("__trop_ctrl_varindex", 2)

    st_numscalar("__trop_panel_varindex", 3)
    st_numscalar("__trop_time_varindex", 4)

    st_numscalar("__trop_n_units", n_units)
    st_numscalar("__trop_n_periods", n_periods)

    st_global("__trop_y_varname", y_varname)
    st_global("__trop_d_varname", d_varname)
    st_global("__trop_panel_varname", panel_varname)
    st_global("__trop_time_varname", time_varname)

    _trop_prepare_time_dist_matrix(n_periods)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_lambda_grids()

  Stores candidate grids for the three regularization parameters
  (lambda_time, lambda_unit, lambda_nn) as Stata matrices, together with
  the corresponding grid lengths as scalars.

  Arguments
    lambda_time_grid   row vector of lambda_time candidates
    lambda_unit_grid   row vector of lambda_unit candidates
    lambda_nn_grid     row vector of lambda_nn   candidates
──────────────────────────────────────────────────────────────────────────────*/

void trop_prepare_lambda_grids(
    real rowvector lambda_time_grid,
    real rowvector lambda_unit_grid,
    real rowvector lambda_nn_grid
)
{
    st_matrix("__trop_lambda_time_grid", lambda_time_grid)
    st_matrix("__trop_lambda_unit_grid", lambda_unit_grid)
    st_matrix("__trop_lambda_nn_grid", lambda_nn_grid)

    st_numscalar("__trop_n_lambda_time", cols(lambda_time_grid))
    st_numscalar("__trop_n_lambda_unit", cols(lambda_unit_grid))
    st_numscalar("__trop_n_lambda_nn", cols(lambda_nn_grid))
}

/*──────────────────────────────────────────────────────────────────────────────
  _trop_prepare_time_dist_matrix()   [private]

  Constructs a T x T symmetric matrix of absolute period distances:

      time_dist[t, s] = |t - s|

  These distances enter the exponential time-decay kernel

      theta_s^{i,t}(lambda) = exp(-lambda_time * |t - s|)

  which down-weights periods far from the target period t.

  Argument
    n_periods   number of time periods (T)
──────────────────────────────────────────────────────────────────────────────*/

void _trop_prepare_time_dist_matrix(real scalar n_periods)
{
    real matrix time_dist
    real scalar t, s

    time_dist = J(n_periods, n_periods, 0)

    for (t = 1; t <= n_periods; t++) {
        for (s = 1; s <= n_periods; s++) {
            time_dist[t, s] = abs(t - s)
        }
    }

    st_matrix("__trop_time_dist", time_dist)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_output_matrices()

  Pre-allocates Stata matrices that the plugin fills via SF_mat_store().
  Must be called before `plugin call`.

  Arguments
    n_units     number of units   (N)
    n_periods   number of periods (T)
    n_treated   number of treated observations (optional; default 1)

  Pre-allocated matrices
    __trop_alpha           N x 1          unit fixed effects   (alpha_i)
    __trop_beta            T x 1          time fixed effects   (beta_t)
    __trop_tau             n_treated x 1  treatment effects    (tau_{i,t})
    __trop_factor_matrix   T x N          low-rank component   (mu)
    __trop_theta           T x 1          time weights   (two-step)
    __trop_omega           N x 1          unit weights   (two-step)
    __trop_delta_time      T x 1          time weights   (joint)
    __trop_delta_unit      N x 1          unit weights   (joint)
──────────────────────────────────────────────────────────────────────────────*/

void trop_prepare_output_matrices(
    real scalar n_units,
    real scalar n_periods,
    | real scalar n_treated
)
{
    if (args() < 3 || n_treated < 1) {
        n_treated = 1
    }

    /* Parameter estimates */
    st_matrix("__trop_alpha", J(n_units, 1, 0))
    st_matrix("__trop_beta", J(n_periods, 1, 0))
    st_matrix("__trop_tau", J(n_treated, 1, 0))
    st_matrix("__trop_factor_matrix", J(n_periods, n_units, 0))

    /* Per-observation diagnostics for the twostep method.  Pre-allocate so
       the plugin can always st_store into these matrices.  -1 sentinels
       distinguish "never written" from "solver returned failure". */
    st_matrix("__trop_converged_by_obs", J(n_treated, 1, -1))
    st_matrix("__trop_n_iters_by_obs",   J(n_treated, 1, -1))

    /* Weight vectors for both estimator variants.
       SF_mat_store() requires the target matrix to exist before the
       plugin call. */
    st_matrix("__trop_theta", J(n_periods, 1, 0))
    st_matrix("__trop_omega", J(n_units, 1, 0))
    st_matrix("__trop_delta_time", J(n_periods, 1, 0))
    st_matrix("__trop_delta_unit", J(n_units, 1, 0))

    /* Pre-allocate gamma output (covariates) as 1 x p row vector.
       trop_prepare_covariates() has already set __trop_n_covariates when
       covariates are present; honour that value.  Fall back to 1 x 1 when
       no covariates so that the matrix always exists for the ADO-side
       `capture confirm matrix __trop_gamma` check.

       Note: st_numscalar() returns J(0,0,.) when the scalar does not exist
       (e.g. after trop_cleanup_temp_vars).  We must read into a matrix
       first to avoid a conformability error assigning void to real scalar. */
    {
        real scalar _p_gamma
        real matrix _tmp_ncov
        _tmp_ncov = st_numscalar("__trop_n_covariates")
        _p_gamma = (rows(_tmp_ncov) > 0 ? _tmp_ncov : .)
        if (_p_gamma >= . | _p_gamma < 1) _p_gamma = 1
        st_matrix("__trop_gamma", J(1, _p_gamma, 0))
    }
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_options()

  Stores estimation options as Stata scalars for the plugin.

  Arguments
    max_iter            maximum iterations for alternating minimization
    tol                 convergence tolerance
    seed                PRNG seed for bootstrap
    bootstrap_reps      number of bootstrap replications
    alpha_level         significance level (e.g. 0.05 for 95% CI)
    verbose             verbosity level (0 = silent)
──────────────────────────────────────────────────────────────────────────────*/

void trop_prepare_options(
    real scalar max_iter,
    real scalar tol,
    real scalar seed,
    real scalar bootstrap_reps,
    real scalar alpha_level,
    real scalar verbose
)
{
    st_numscalar("__trop_max_iter", max_iter)
    st_numscalar("__trop_tol", tol)
    st_numscalar("__trop_seed", seed)
    st_numscalar("__trop_n_bootstrap", bootstrap_reps)
    st_numscalar("__trop_alpha_level", alpha_level)
    st_numscalar("__trop_verbose", verbose)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_bootstrap()

  Stores bootstrap parameters and pre-allocates the output matrix that
  collects replicate estimates.

  Arguments
    n_bootstrap   number of bootstrap replications
    alpha         significance level for confidence intervals
    seed          PRNG seed
    lambda_time   selected lambda_time
    lambda_unit   selected lambda_unit
    lambda_nn     selected lambda_nn
    max_iter      maximum iterations per replicate
    tol           convergence tolerance per replicate

  Optional argument
    ddof          variance denominator selector forwarded to the plugin:
                  1 = sample variance 1/(B-1) (default);
                  0 = paper Algorithm 3 population variance 1/B.
                  Any other value collapses to 1.  When omitted the scalar
                  __trop_bs_ddof is left unset and the plugin applies its
                  own default (1), preserving pre-existing call sites.

  Pre-allocated matrix
    __trop_bootstrap_estimates   n_bootstrap x 1, initialised to missing
──────────────────────────────────────────────────────────────────────────────*/

void trop_prepare_bootstrap(
    real scalar n_bootstrap,
    real scalar alpha,
    real scalar seed,
    real scalar lambda_time,
    real scalar lambda_unit,
    real scalar lambda_nn,
    real scalar max_iter,
    real scalar tol,
    | real scalar ddof
)
{
    real scalar ddof_eff

    st_numscalar("__trop_n_bootstrap", n_bootstrap)
    /* Named __trop_bs_alpha to avoid collision with the unit-fixed-effects
       matrix __trop_alpha. */
    st_numscalar("__trop_bs_alpha", alpha)
    st_numscalar("__trop_seed", seed)

    st_numscalar("__trop_lambda_time", lambda_time)
    st_numscalar("__trop_lambda_unit", lambda_unit)
    st_numscalar("__trop_lambda_nn", lambda_nn)

    st_numscalar("__trop_max_iter", max_iter)
    st_numscalar("__trop_tol", tol)

    if (args() >= 9 & ddof < .) {
        /* Clamp to {0, 1}; any other value collapses to 1. */
        ddof_eff = (ddof == 0) ? 0 : 1
        st_numscalar("__trop_bs_ddof", ddof_eff)
    }

    st_matrix("__trop_bootstrap_estimates", J(n_bootstrap, 1, .))
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_covariates()

  Reads covariate variables from the Stata dataset, pivots each into a
  (T x N) matrix (same layout as Y), then stacks them column-wise into
  a (T*N) x p matrix stored as __trop_covariates.  The row ordering is
  time-major: row k = t * n_units + i, matching the Rust internal layout.

  Arguments
    cov_varnames    row vector of covariate variable names
    panel_idx_var   panel index variable (1..N from egen group)
    time_idx_var    time index variable (1..T from egen group)
    touse_var       estimation-sample marker
    n_units         number of cross-sectional units (N)
    n_periods       number of time periods (T)

  Stored matrices
    __trop_covariates   (T*N) x p covariate matrix

  Stored scalars
    __trop_n_covariates   number of covariates (p)

  Returns
    0 on success; nonzero on failure
──────────────────────────────────────────────────────────────────────────────*/

real scalar trop_prepare_covariates(
    string rowvector cov_varnames,
    string scalar panel_idx_var,
    string scalar time_idx_var,
    string scalar touse_var,
    real scalar n_units,
    real scalar n_periods
)
{
    real scalar p, n_obs, i, obs, t_val, u_val, row_idx
    real colvector panel_idx, time_idx, cov_data
    real matrix X
    string scalar varname

    p = cols(cov_varnames)
    n_obs = n_periods * n_units

    /* Read panel and time indices */
    panel_idx = st_data(., panel_idx_var, touse_var)
    time_idx  = st_data(., time_idx_var,  touse_var)

    /* Allocate output matrix */
    X = J(n_obs, p, 0)

    /* Fill each covariate column */
    for (i = 1; i <= p; i++) {
        varname = cov_varnames[i]
        cov_data = st_data(., varname, touse_var)

        /* Map observation-level data to (t-1)*n_units + (u-1) + 1 row */
        for (obs = 1; obs <= rows(cov_data); obs++) {
            t_val = time_idx[obs]    /* 1-based time index */
            u_val = panel_idx[obs]   /* 1-based unit index */
            /* Row index: time-major, 1-based for Mata */
            row_idx = (t_val - 1) * n_units + u_val
            if (row_idx >= 1 & row_idx <= n_obs) {
                X[row_idx, i] = cov_data[obs]
            }
        }
    }

    /* Secondary validation: reject if covariate matrix contains missing values.
       This should never trigger because the ADO layer rejects missing covariates
       upstream; included as defensive programming. */
    if (hasmissing(X)) {
        _error(416, "internal error: covariate matrix contains missing values after validation")
    }

    /* Store to Stata */
    st_matrix("__trop_covariates", X)
    st_numscalar("__trop_n_covariates", p)

    /* Pre-allocate gamma output as 1 x p row vector (Stata convention).
       The plugin writes SF_mat_store("__trop_gamma", 1, j+1, val) so this
       must exist with the correct dimensions before the plugin call. */
    st_matrix("__trop_gamma", J(1, p, 0))

    return(0)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_pweights()

  Extracts a strictly positive, constant-within-unit pweight vector from a
  Stata variable and stores the resulting N x 1 column of per-unit weights
  in the matrix __trop_unit_weights.  Also sets __trop_use_weights = 1 so
  the plugin dispatches to the weighted Rust ABI.

  Arguments
    weight_var       name of the pweight variable
    panel_idx_var    unit identifier (1..N) after egen group
    touse_var        estimation-sample marker
    n_units          number of units (N)

  Returns
    0 on success; nonzero Stata return code on validation failure:
      198 if any weight is missing, non-positive, or non-finite
      198 if pweight is not constant within a unit

  Side effects
    __trop_unit_weights  (N x 1 matrix)
    __trop_use_weights   = 1
──────────────────────────────────────────────────────────────────────────────*/

real scalar trop_prepare_pweights(
    string scalar weight_var,
    string scalar panel_idx_var,
    string scalar touse_var,
    real scalar n_units
)
{
    real colvector panel_idx, w_obs, w_unit, seen
    real scalar i, idx, wi

    panel_idx = st_data(., panel_idx_var, touse_var)
    w_obs     = st_data(., weight_var,    touse_var)

    if (rows(panel_idx) != rows(w_obs)) {
        errprintf("pweight (%s) and panel index lengths differ\n", weight_var)
        return(459)
    }

    /* Any missing / non-positive pweight cell is a hard error. */
    for (i = 1; i <= rows(w_obs); i++) {
        wi = w_obs[i]
        if (wi >= . | wi <= 0) {
            errprintf("pweight %s must be strictly positive; found %g at obs %g\n",
                      weight_var, wi, i)
            return(459)
        }
    }

    /* Collect the first-seen weight per unit, then enforce that every
       subsequent observation in the unit reports the same value. */
    w_unit = J(n_units, 1, .)
    seen   = J(n_units, 1, 0)

    for (i = 1; i <= rows(w_obs); i++) {
        idx = panel_idx[i]
        if (idx < 1 | idx > n_units) {
            errprintf("panel index %g out of range [1, %g]\n", idx, n_units)
            return(459)
        }
        if (!seen[idx]) {
            w_unit[idx] = w_obs[i]
            seen[idx]   = 1
        }
        else if (reldif(w_unit[idx], w_obs[i]) > 1e-12) {
            errprintf("pweight %s is not constant within unit %g (found %g and %g)\n",
                      weight_var, idx, w_unit[idx], w_obs[i])
            return(459)
        }
    }

    /* Units absent from the touse sample receive 0; the Rust aggregator
       ignores non-positive weights, so this is safe. */
    for (i = 1; i <= n_units; i++) {
        if (!seen[i]) w_unit[i] = 0
    }

    st_matrix("__trop_unit_weights", w_unit)
    st_numscalar("__trop_use_weights", 1)
    st_global("__trop_weight_var", weight_var)

    return(0)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_prepare_survey_design()

  Reads survey design variables (strata, PSU, FPC) from the Stata dataset,
  collapses observation-level values to unit-level (verifying within-unit
  constancy), encodes them as 0-based integer indices, and stores the
  results as Stata matrices for the Rao-Wu bootstrap plugin.

  Arguments
    strata_var      name of the stratification variable
    psu_var         name of the primary sampling unit variable
    fpc_var         name of the FPC variable ("" if not specified)
    panel_idx_var   panel index variable (1..N, from egen group)
    touse_var       estimation-sample marker
    n_units         number of cross-sectional units (N)
    do_nest         1 = nest PSU within strata (combine strata*PSU labels)

  Stored matrices
    __trop_strata   N x 1  integer-encoded strata (0-based)
    __trop_psu      N x 1  integer-encoded PSU    (0-based)
    __trop_fpc      N x 1  FPC values (or empty if fpc_var=="")

  Stored scalars
    __trop_has_survey_design = 1
    __trop_n_strata          number of distinct strata
    __trop_n_psu             number of distinct PSU

  Returns
    0 on success; nonzero Stata return code on validation failure.
──────────────────────────────────────────────────────────────────────────────*/

real scalar trop_prepare_survey_design(
    string scalar strata_var,
    string scalar psu_var,
    string scalar fpc_var,
    string scalar panel_idx_var,
    string scalar touse_var,
    real scalar n_units,
    real scalar do_nest
)
{
    real colvector panel_idx, strata_obs, psu_obs, fpc_obs
    real colvector strata_unit, psu_unit, fpc_unit, seen
    real colvector strata_encoded, psu_encoded
    real scalar i, idx, n_obs
    real scalar n_strata, n_psu
    real colvector strata_labels, psu_labels
    real scalar found

    /* 1. Read observation-level data */
    panel_idx  = st_data(., panel_idx_var, touse_var)
    strata_obs = st_data(., strata_var,    touse_var)
    psu_obs    = st_data(., psu_var,       touse_var)
    n_obs = rows(panel_idx)

    if (fpc_var != "") {
        fpc_obs = st_data(., fpc_var, touse_var)
    }

    /* 2. Collapse to unit-level and verify within-unit constancy */
    strata_unit = J(n_units, 1, .)
    psu_unit    = J(n_units, 1, .)
    fpc_unit    = J(n_units, 1, .)
    seen        = J(n_units, 1, 0)

    for (i = 1; i <= n_obs; i++) {
        idx = panel_idx[i]
        if (idx < 1 | idx > n_units) {
            errprintf("panel index %g out of range [1, %g]\n", idx, n_units)
            return(459)
        }
        if (!seen[idx]) {
            strata_unit[idx] = strata_obs[i]
            psu_unit[idx]    = psu_obs[i]
            if (fpc_var != "") fpc_unit[idx] = fpc_obs[i]
            seen[idx] = 1
        }
        else {
            /* Verify within-unit constancy */
            if (strata_unit[idx] != strata_obs[i]) {
                errprintf("strata variable '%s' is not constant within unit %g "
                    + "(found %g and %g)\n",
                    strata_var, idx, strata_unit[idx], strata_obs[i])
                return(459)
            }
            if (psu_unit[idx] != psu_obs[i]) {
                errprintf("psu variable '%s' is not constant within unit %g "
                    + "(found %g and %g)\n",
                    psu_var, idx, psu_unit[idx], psu_obs[i])
                return(459)
            }
            if (fpc_var != "") {
                if (reldif(fpc_unit[idx], fpc_obs[i]) > 1e-12) {
                    errprintf("fpc variable '%s' is not constant within unit %g "
                        + "(found %g and %g)\n",
                        fpc_var, idx, fpc_unit[idx], fpc_obs[i])
                    return(459)
                }
            }
        }
    }

    /* 3. Integer encoding (0-based) */
    /* If nest=1, combine strata+PSU into a single label before encoding PSU */
    if (do_nest) {
        /* Create composite label: strata * 1e9 + psu */
        for (i = 1; i <= n_units; i++) {
            psu_unit[i] = strata_unit[i] * 1e9 + psu_unit[i]
        }
    }

    /* Encode strata: map unique values to 0, 1, 2, ... */
    strata_labels = J(0, 1, .)
    strata_encoded = J(n_units, 1, 0)
    for (i = 1; i <= n_units; i++) {
        found = 0
        for (idx = 1; idx <= rows(strata_labels); idx++) {
            if (strata_labels[idx] == strata_unit[i]) {
                strata_encoded[i] = idx - 1
                found = 1
                break
            }
        }
        if (!found) {
            strata_labels = strata_labels \ strata_unit[i]
            strata_encoded[i] = rows(strata_labels) - 1
        }
    }
    n_strata = rows(strata_labels)

    /* Encode PSU: map unique values to 0, 1, 2, ... */
    psu_labels = J(0, 1, .)
    psu_encoded = J(n_units, 1, 0)
    for (i = 1; i <= n_units; i++) {
        found = 0
        for (idx = 1; idx <= rows(psu_labels); idx++) {
            if (psu_labels[idx] == psu_unit[i]) {
                psu_encoded[i] = idx - 1
                found = 1
                break
            }
        }
        if (!found) {
            psu_labels = psu_labels \ psu_unit[i]
            psu_encoded[i] = rows(psu_labels) - 1
        }
    }
    n_psu = rows(psu_labels)

    /* 3a. Validate FPC >= n_psu_in_stratum for each stratum */
    if (fpc_var != "") {
        real scalar h, n_psu_h, fpc_h
        real colvector psu_in_h
        for (h = 0; h < n_strata; h++) {
            psu_in_h = J(0, 1, .)
            fpc_h = .
            for (i = 1; i <= n_units; i++) {
                if (strata_encoded[i] == h) {
                    if (fpc_h == .) fpc_h = fpc_unit[i]
                    found = 0
                    for (idx = 1; idx <= rows(psu_in_h); idx++) {
                        if (psu_in_h[idx] == psu_encoded[i]) {
                            found = 1
                            break
                        }
                    }
                    if (!found) {
                        psu_in_h = psu_in_h \ psu_encoded[i]
                    }
                }
            }
            n_psu_h = rows(psu_in_h)
            if (fpc_h < n_psu_h) {
                errprintf("fpc() value (%.0f) is less than the number of "
                    + "PSUs (%g) in stratum %g\n", fpc_h, n_psu_h, h + 1)
                return(459)
            }
        }
    }

    /* 3b. nest=0: verify PSU codes are globally unique across strata */
    if (!do_nest) {
        real colvector psu_unique, strata_for_p, distinct_s
        real scalar p_idx, p_val
        psu_unique = uniqrows(psu_encoded)
        for (p_idx = 1; p_idx <= rows(psu_unique); p_idx++) {
            p_val = psu_unique[p_idx]
            strata_for_p = select(strata_encoded, psu_encoded :== p_val)
            distinct_s = uniqrows(strata_for_p)
            if (rows(distinct_s) > 1) {
                errprintf("psu variable is not globally unique "
                    + "(nest option not specified): ")
                errprintf("PSU code %g appears in %g different strata\n",
                    p_val, rows(distinct_s))
                errprintf("  Use the 'nest' option if PSUs are nested "
                    + "within strata\n")
                return(459)
            }
        }
    }

    /* 4. Store to Stata matrices */
    st_matrix("__trop_strata", strata_encoded)
    st_matrix("__trop_psu", psu_encoded)
    if (fpc_var != "") {
        st_matrix("__trop_fpc", fpc_unit)
    }
    else {
        /* Empty 0x1 sentinel: plugin checks rows==0 to skip FPC */
        st_matrix("__trop_fpc", J(0, 1, .))
    }

    /* 5. Store metadata scalars */
    st_numscalar("__trop_has_survey_design", 1)
    st_numscalar("__trop_n_strata", n_strata)
    st_numscalar("__trop_n_psu", n_psu)

    return(0)
}

/*──────────────────────────────────────────────────────────────────────────────
  trop_cleanup_temp_vars()

  Drops every __trop_* scalar, matrix, and global macro from the Stata
  namespace.  Uses st_dir() to enumerate names programmatically so that
  newly introduced temporaries are covered without code changes.
──────────────────────────────────────────────────────────────────────────────*/

void trop_cleanup_temp_vars()
{
    string colvector names
    real scalar i

    /* Scalars */
    names = st_dir("global", "numscalar", "__trop_*")
    for (i = 1; i <= length(names); i++) {
        stata("capture scalar drop " + names[i])
    }

    /* Matrices */
    names = st_dir("global", "matrix", "__trop_*")
    for (i = 1; i <= length(names); i++) {
        stata("capture matrix drop " + names[i])
    }

    /* Global macros */
    names = st_dir("global", "macro", "__trop_*")
    for (i = 1; i <= length(names); i++) {
        st_global(names[i], "")
    }
    st_global("__trop_touse_var", "")
    st_global("__trop_varlist", "")
    st_global("__trop_y_varname", "")
    st_global("__trop_d_varname", "")
    st_global("__trop_panel_varname", "")
    st_global("__trop_time_varname", "")

    /* Clean up the shared random seed scalar (not __trop_ prefixed) */
    stata("capture scalar drop TROP_GLOBAL_SEED")
}

end