libqe coverage - 71.20%

Files
Source

/**
 * @file adjacency.hpp
 * @brief Combinatorics and vector/matrix utilities for ENA connection networks.
 *
 * Pure C++ / Armadillo — no Rcpp dependency.
 * Include @c <RcppArmadillo.h> before this header when building inside an R package.
 */
#ifndef LIBQE_ADJACENCY_HPP
#define LIBQE_ADJACENCY_HPP

#include <armadillo>
#include <string>
#include <vector>
#include <cmath>

namespace qe {

// ---------------------------------------------------------------------------
/// @name Combinatorics
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Number of unordered pairs from @p n items (n choose 2).
 *
 * @param n Number of codes/nodes.
 * @returns @c (n * (n - 1)) / 2.
 */
inline int choose_two(int n) {
    return (n * (n - 1)) / 2;
}

/**
 * @brief Number of undirected connection columns for @p n codes.
 *
 * @param n        Number of codes/nodes.
 * @param diagonal Include self-connections (the diagonal, @c A & A).
 * @returns @c choose_two(n), or @c n(n+1)/2 with @p diagonal.
 */
inline int n_connections(int n, bool diagonal = false) {
    return diagonal ? (n * (n + 1)) / 2 : choose_two(n);
}

/**
 * @brief Upper-triangle index pairs for a symmetric matrix of side @p len.
 *
 * Column-major order: (0,1), (0,2), (1,2), ... -- or, with @p diagonal,
 * (0,0), (0,1), (1,1), (0,2), (1,2), (2,2), ... (Armadillo's
 * @c trimatu_ind order), so a self-connection's row and column index are
 * equal.
 *
 * @param len      Side length of the square matrix (number of codes).
 * @param row      Controls what is returned:
 *   - @c -1 (default): 2 × k matrix of [row_idx; col_idx].
 *   - @c  0: row indices only (1 × k).
 *   - @c  1: column indices only (1 × k).
 * @param diagonal Include the diagonal (self-connections).  Default @c false.
 * @returns An @c arma::umat of shape 2×k or 1×k depending on @p row, with
 *          @c k = n_connections(len, diagonal).
 * @note Equivalent to @c triIndices() in rENA/ena.cpp and tma/code.cpp
 *       (without @p diagonal).
 */
inline arma::umat connection_indices(int len, int row = -1, bool diagonal = false) {
    int vS = n_connections(len, diagonal);
    int s  = 0;

    arma::umat vR(2, vS, arma::fill::zeros);
    arma::umat vRone(1, vS, arma::fill::zeros);

    for (int i = 1; i <= len; i++) {
        for (int j = 0; j < (diagonal ? i : i - 1); j++) {
            vR(0, s) = j;
            vR(1, s) = i - 1;
            if (row == 0)      vRone[s] = j;
            else if (row == 1) vRone[s] = i - 1;
            s++;
        }
    }
    return (row == -1) ? vR : vRone;
}

/// @}

// ---------------------------------------------------------------------------
/// @name Vector ↔ upper-triangle conversions
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Compute pairwise products and return as a flat upper-triangle vector.
 *
 * For each pair @c (j, i) with @c j < i, outputs @c v[j] * v[i].
 * The result has length @c choose_two(v.size()).
 *
 * @param v A 1 × p row vector (or p-element matrix).
 * @returns A row vector of length @c choose_two(p).
 * @note Equivalent to @c vector_to_ut() in rENA/ena.cpp.
 */
inline arma::rowvec code_connections(arma::mat v) {
    int vL = v.size();
    int vS = choose_two(vL);
    arma::rowvec out(vS, arma::fill::zeros);
    int s = 0;
    for (int i = 2; i <= vL; i++) {
        for (int j = 0; j < i - 1; j++) {
            out[s] = v[j] * v[i - 1];
            s++;
        }
    }
    return out;
}

/**
 * @brief Fold a directed n² vector into an undirected upper-triangle vector.
 *
 * Symmetric elements are summed: A→B + B→A for every pair.  With
 * @p diagonal the self-connections (A→A) are kept, each once, in the
 * @c connection_indices(n, -1, true) order (A & A, A & B, B & B, ...).
 *
 * @param v        Flat directed vector of length n², column-major
 *                 (n = sqrt(v.size())).
 * @param diagonal Keep the diagonal.  Default @c false.
 * @returns Upper-triangle row vector of length @c n_connections(n, diagonal).
 * @note Equivalent to @c vector_to_summed_uppertri() in tma/code.cpp
 *       (without @p diagonal).
 */
inline arma::rowvec fold_directed_network(arma::vec v, bool diagonal = false) {
    int n    = static_cast<int>(std::round(std::sqrt(static_cast<double>(v.size()))));

    arma::mat m(v);
    m.reshape(n, n);
    arma::mat folded = m + m.t();
    if (diagonal) folded.diag() = m.diag();   // a self-connection is not a pair

    arma::colvec flat     = folded.as_col();
    arma::uvec   up_inds  = arma::trimatu_ind(arma::size(folded), diagonal ? 0 : 1);
    return arma::conv_to<arma::rowvec>::from(flat.elem(up_inds));
}

/**
 * @brief Flatten an adjacency matrix to a vector.
 *
 * @param x        Square adjacency matrix (n × n).
 * @param full     If @c true (default), returns the full n² vector (directed).
 *                 If @c false, returns the upper-triangle only (undirected).
 * @param diagonal With @p full @c false: include the diagonal, in the
 *                 @c connection_indices(n, -1, true) order.  Default @c false.
 * @returns Row vector of length n² or @c n_connections(n, diagonal).
 * @note Equivalent to @c adjacency_matrix_to_vector() in tma/code.cpp
 *       (without @p diagonal).
 */
inline arma::rowvec network_to_vector(arma::mat x, bool full = true, bool diagonal = false) {
    if (full) return arma::vectorise(x).t();
    if (diagonal) {
        arma::uvec up_inds = arma::trimatu_ind(arma::size(x), 0);
        return arma::vectorise(x(up_inds)).t();
    }
    arma::mat combined  = arma::trimatu(x, 1) + arma::trimatl(x, -1);
    arma::uvec up_inds  = arma::trimatu_ind(arma::size(combined), 1);
    return arma::vectorise(combined(up_inds)).t();
}

/// @}

// ---------------------------------------------------------------------------
/// @name String utilities
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Generate "A & B" pair names for every upper-triangle position.
 *
 * For code names @c {"A", "B", "C"} returns @c {"A & B", "A & C", "B & C"};
 * with @p diagonal, @c {"A & A", "A & B", "B & B", "A & C", "B & C", "C & C"}
 * (the @c connection_indices(p, -1, true) order).
 *
 * @param v        Ordered list of code names (length p).
 * @param diagonal Include self-connections.  Default @c false.
 * @returns Vector of @c n_connections(p, diagonal) label strings.
 * @note Equivalent to @c svector_to_ut() in rENA/ena.cpp (without @p diagonal).
 */
inline std::vector<std::string> connection_names(std::vector<std::string> v,
                                                 bool diagonal = false) {
    int vL = v.size();
    int vS = n_connections(vL, diagonal);
    std::vector<std::string> out(vS);
    int s = 0;
    for (int i = 1; i <= vL; i++) {
        for (int j = 0; j < (diagonal ? i : i - 1); j++) {
            out[s] = v[j] + " & " + v[i - 1];
            s++;
        }
    }
    return out;
}

/// @}

} // namespace qe

#endif // LIBQE_ADJACENCY_HPP

/** @file door.hpp
 *  @brief Door temporal pooling and smoothing functions for trajectory modeling.
 *
 *  Provides sliding lookback window pooling (with uniform or linear weighting,
 *  sum or mean aggregation, and optional segment boundary resets) as well as
 *  Exponential Moving Average (EMA) smoothing over consecutive time steps.
 */
#ifndef LIBQE_DOOR_HPP
#define LIBQE_DOOR_HPP

#include <armadillo>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>

namespace qe {

/** @brief Apply a sliding lookback window pool over a single unit's connection matrix.
 *
 *  For each row @c i, pools rows in the window [@c max(0, i - lookback_size + 1), @c i].
 *  If @p segment_ids is non-empty, the window does not cross rows with a different segment ID.
 *
 *  @param[in] block            Matrix of connection counts for one unit (n_rows × n_features).
 *  @param[in] lookback_size    Number of rows in the lookback window (>= 1).
 *  @param[in] aggregate_mean   When @c true, averages rows; when @c false, sums rows.
 *  @param[in] weighting_linear When @c true, applies linearly increasing weights
 *                              (oldest row in window has weight 1, current row has weight k).
 *                              When @c false, all rows in window have equal weight 1.
 *  @param[in] segment_ids      Optional vector of segment/activity IDs (length n_rows).
 *                              Window pooling resets at segment transitions.
 *
 *  @returns A matrix of shape n_rows × n_features containing door-smoothed values.
 */
inline arma::mat lookback_block(
    const arma::mat& block,
    int lookback_size = 20,
    bool aggregate_mean = false,
    bool weighting_linear = false,
    const std::vector<int>& segment_ids = {}
) {
    int n_rows = static_cast<int>(block.n_rows);
    int n_cols = static_cast<int>(block.n_cols);
    arma::mat out(n_rows, n_cols, arma::fill::zeros);

    if (n_rows == 0 || n_cols == 0) return out;
    int k_max = std::max(1, lookback_size);
    bool has_segments = (!segment_ids.empty() && static_cast<int>(segment_ids.size()) >= n_rows);

    for (int i = 0; i < n_rows; ++i) {
        int start_idx = std::max(0, i - k_max + 1);

        if (has_segments) {
            int cur_seg = segment_ids[i];
            for (int s = i; s >= start_idx; --s) {
                if (segment_ids[s] != cur_seg) {
                    start_idx = s + 1;
                    break;
                }
            }
        }

        int win_len = i - start_idx + 1;
        if (win_len <= 0) continue;

        arma::rowvec acc(n_cols, arma::fill::zeros);
        arma::rowvec denom(n_cols, arma::fill::zeros);

        for (int j = start_idx; j <= i; ++j) {
            double w = 1.0;
            if (weighting_linear) {
                // Linear weight: 1 for earliest in window, win_len for current row i
                w = static_cast<double>(j - start_idx + 1);
            }
            for (int col = 0; col < n_cols; ++col) {
                double value = block(j, col);
                if (!std::isnan(value)) {
                    acc[col] += w * value;
                    denom[col] += w;
                }
            }
        }

        if (aggregate_mean) {
            for (int col = 0; col < n_cols; ++col) {
                acc[col] = (denom[col] > 0.0) ? (acc[col] / denom[col]) : arma::datum::nan;
            }
        }

        out.row(i) = acc;
    }

    return out;
}

/** @brief Apply Exponential Moving Average (EMA) smoothing over a single unit's connection matrix.
 *
 *  Computes S_0 = Y_0, S_t = alpha * Y_t + (1 - alpha) * S_{t-1}.
 *  If @p segment_ids is non-empty, S_t resets to Y_t when transitioning to a new segment.
 *
 *  @param[in] block       Matrix of connection counts for one unit (n_rows × n_features).
 *  @param[in] alpha       Smoothing parameter alpha in (0, 1].
 *  @param[in] segment_ids Optional vector of segment/activity IDs (length n_rows).
 *
 *  @returns A matrix of shape n_rows × n_features containing EMA-smoothed values.
 */
inline arma::mat ema_block(
    const arma::mat& block,
    double alpha = 0.1,
    const std::vector<int>& segment_ids = {}
) {
    int n_rows = static_cast<int>(block.n_rows);
    int n_cols = static_cast<int>(block.n_cols);
    arma::mat out(n_rows, n_cols, arma::fill::zeros);

    if (n_rows == 0 || n_cols == 0) return out;
    double a = std::max(1e-6, std::min(1.0, alpha));
    bool has_segments = (!segment_ids.empty() && static_cast<int>(segment_ids.size()) >= n_rows);

    out.row(0) = block.row(0);

    for (int i = 1; i < n_rows; ++i) {
        bool reset = has_segments && (segment_ids[i] != segment_ids[i - 1]);
        if (reset) {
            out.row(i) = block.row(i);
        } else {
            for (int col = 0; col < n_cols; ++col) {
                double cur = block(i, col);
                if (std::isnan(cur)) {
                    cur = out(i - 1, col);
                }
                out(i, col) = a * cur + (1.0 - a) * out(i - 1, col);
            }
        }
    }

    return out;
}

/** @brief Apply lookback door pooling to an entire connection matrix across multiple units.
 *
 *  @param[in] conn_counts        Complete connection matrix (all units).
 *  @param[in] unit_row_indices   Vector of row-index vectors, one per unit (0-indexed).
 *  @param[in] lookback_sizes     Vector of lookback sizes, one per unit.
 *  @param[in] aggregate_mean     When @c true, compute mean; when @c false, compute sum.
 *  @param[in] weighting_linear   When @c true, use linear weighting.
 *  @param[in] segment_ids        Optional vector of segment IDs for every row in conn_counts.
 *
 *  @returns Smoothed connection matrix with identical dimensions to @p conn_counts.
 */
inline arma::mat apply_door_lookback(
    const arma::mat& conn_counts,
    const std::vector<std::vector<int>>& unit_row_indices,
    const std::vector<int>& lookback_sizes,
    bool aggregate_mean = false,
    bool weighting_linear = false,
    const std::vector<int>& segment_ids = {}
) {
    arma::mat out = conn_counts;
    size_t n_units = unit_row_indices.size();
    bool has_segments = (!segment_ids.empty() && segment_ids.size() >= conn_counts.n_rows);

    for (size_t u = 0; u < n_units; ++u) {
        const auto& rows = unit_row_indices[u];
        if (rows.empty()) continue;

        int k = (u < lookback_sizes.size()) ? lookback_sizes[u] : 20;
        int n_u = static_cast<int>(rows.size());

        arma::mat block(n_u, conn_counts.n_cols);
        std::vector<int> block_segs(n_u, 0);

        for (int i = 0; i < n_u; ++i) {
            int r = rows[i];
            block.row(i) = conn_counts.row(r);
            if (has_segments) block_segs[i] = segment_ids[r];
        }

        arma::mat smoothed = lookback_block(block, k, aggregate_mean, weighting_linear,
                                            has_segments ? block_segs : std::vector<int>{});

        for (int i = 0; i < n_u; ++i) {
            out.row(rows[i]) = smoothed.row(i);
        }
    }

    return out;
}

/** @brief Apply EMA door smoothing to an entire connection matrix across multiple units.
 */
inline arma::mat apply_door_ema(
    const arma::mat& conn_counts,
    const std::vector<std::vector<int>>& unit_row_indices,
    double alpha = 0.1,
    const std::vector<int>& segment_ids = {}
) {
    arma::mat out = conn_counts;
    size_t n_units = unit_row_indices.size();
    bool has_segments = (!segment_ids.empty() && segment_ids.size() >= conn_counts.n_rows);

    for (size_t u = 0; u < n_units; ++u) {
        const auto& rows = unit_row_indices[u];
        if (rows.empty()) continue;

        int n_u = static_cast<int>(rows.size());
        arma::mat block(n_u, conn_counts.n_cols);
        std::vector<int> block_segs(n_u, 0);

        for (int i = 0; i < n_u; ++i) {
            int r = rows[i];
            block.row(i) = conn_counts.row(r);
            if (has_segments) block_segs[i] = segment_ids[r];
        }

        arma::mat smoothed = ema_block(block, alpha, has_segments ? block_segs : std::vector<int>{});

        for (int i = 0; i < n_u; ++i) {
            out.row(rows[i]) = smoothed.row(i);
        }
    }

    return out;
}

} // namespace qe

#endif // LIBQE_DOOR_HPP

/** @file linalg_fallback.hpp
 *  @brief Thin LAPACK-free implementations of SVD, QR, and symmetric
 *         eigendecomposition for use in WASM / no-LAPACK builds.
 *
 *  When @c ARMA_USE_LAPACK is defined these wrappers delegate directly to the
 *  standard Armadillo routines (@c arma::svd, @c arma::qr, @c arma::eig_sym),
 *  so there is zero overhead in normal builds.  When LAPACK is not available
 *  (e.g. Emscripten WASM builds) the wrappers fall back to pure C++
 *  implementations that work for the small matrices typical in ENA
 *  (p ≤ ~100 connection dimensions).
 *
 *  All routines live in namespace @c qe::linalg and operate on @c arma::mat
 *  (double precision).
 *
 *  Fallback algorithms:
 *  - **SVD**      — one-sided Jacobi iterations on \f$A^T\f$ (Demmel & Veselić 1992)
 *  - **QR**       — modified Gram-Schmidt
 *  - **eig_sym**  — cyclic Jacobi (classical Jacobi sweep)
 */

#ifndef LIBQE_LINALG_FALLBACK_HPP
#define LIBQE_LINALG_FALLBACK_HPP

#include <armadillo>
#include <cmath>
#include <limits>

namespace qe {
namespace linalg {

/// @cond INTERNAL
// ── helpers ───────────────────────────────────────────────────────────────────

namespace detail {

/// @brief Modified Gram-Schmidt orthonormalization applied in-place.
///
/// Columns of @p Q are successively orthonormalized against all previously
/// processed columns.  A column whose norm falls below 1e-14 after projection
/// is left unchanged (not zeroed) to avoid introducing NaNs.
///
/// @param Q  Matrix whose columns are orthonormalized in-place (modified).
inline void mgs(arma::mat& Q) {
    const arma::uword p = Q.n_cols;
    for (arma::uword j = 0; j < p; ++j) {
        for (arma::uword i = 0; i < j; ++i)
            Q.col(j) -= arma::dot(Q.col(i), Q.col(j)) * Q.col(i);
        double n = arma::norm(Q.col(j));
        if (n > 1e-14) Q.col(j) /= n;
    }
}

/// @brief One cyclic Jacobi sweep on a p×p symmetric matrix.
///
/// Iterates over all super-diagonal index pairs @c (q,r) and applies a
/// two-sided Jacobi rotation that annihilates @c A(q,r).  The same rotation
/// is accumulated in @p V so that, starting with @c V = I, the full sequence
/// of sweeps yields @c V = eigenvector matrix.
///
/// @param A  Symmetric p×p matrix, updated in-place toward diagonal form.
/// @param V  Rotation accumulator (p×p); initialize to identity to recover
///           eigenvectors.
///
/// @note Algorithm follows the classical cyclic-by-rows Jacobi scheme.
///       Off-diagonal entries no larger than 1e-14 × (|a_qq| + |a_rr|) are
///       skipped for numerical stability.  The `<=` (not `<`) guard correctly
///       handles the all-zero case (a_qq = a_rr = a_qr = 0) that would
///       otherwise produce NaN via 0/0 in the theta computation.
inline void jacobi_sweep(arma::mat& A, arma::mat& V) {
    const arma::uword p = A.n_rows;
    for (arma::uword q = 0; q < p - 1; ++q) {
        for (arma::uword r = q + 1; r < p; ++r) {
            double aqq = A(q, q), arr = A(r, r), aqr = A(q, r);
            if (std::abs(aqr) <= 1e-14 * (std::abs(aqq) + std::abs(arr)))
                continue;
            double theta = (arr - aqq) / (2.0 * aqr);
            double t = (theta >= 0.0)
                ? 1.0 / (theta + std::sqrt(1.0 + theta * theta))
                : 1.0 / (theta - std::sqrt(1.0 + theta * theta));
            double c = 1.0 / std::sqrt(1.0 + t * t); ///< cosine of Jacobi rotation angle
            double s = t * c;                          ///< sine of Jacobi rotation angle
            // Apply Jacobi rotation to A (symmetric update)
            double new_qq = aqq - t * aqr;
            double new_rr = arr + t * aqr;
            A(q, q) = new_qq;  A(r, r) = new_rr;  A(q, r) = A(r, q) = 0.0;
            for (arma::uword k = 0; k < p; ++k) {
                if (k == q || k == r) continue;
                double akq = A(k, q), akr = A(k, r);
                A(k, q) = A(q, k) = c * akq - s * akr;
                A(k, r) = A(r, k) = s * akq + c * akr;
            }
            // Accumulate rotation in V
            for (arma::uword k = 0; k < p; ++k) {
                double vkq = V(k, q), vkr = V(k, r);
                V(k, q) = c * vkq - s * vkr;
                V(k, r) = s * vkq + c * vkr;
            }
        }
    }
}

} // namespace detail

/// @endcond

/// @name Eigendecomposition
/// @{

// ── eig_sym ───────────────────────────────────────────────────────────────────

/// @brief Eigendecomposition of a real symmetric matrix.
///
/// Computes eigenvalues and eigenvectors of the symmetric p×p matrix @p S and
/// returns them sorted in **ascending** eigenvalue order.
///
/// @param[out] eigval  p×1 vector of eigenvalues, ascending.
/// @param[out] eigvec  p×p matrix whose columns are the corresponding
///                     orthonormal eigenvectors.
/// @param[in]  S       Real symmetric p×p input matrix (not modified).
///
/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::eig_sym
///          directly.  Without LAPACK the fallback runs up to @c 50*p cyclic
///          Jacobi sweeps; convergence is guaranteed for well-conditioned
///          matrices but may be slow for p > 200.
///
/// @note Fallback algorithm: classical cyclic-by-rows Jacobi iteration.
///       Convergence is declared when the Frobenius norm of the strictly
///       off-diagonal part drops below 1e-13.
inline void eig_sym(arma::vec& eigval, arma::mat& eigvec, const arma::mat& S) {
#if defined(ARMA_USE_LAPACK)
    arma::eig_sym(eigval, eigvec, S);
#else
    const arma::uword p = S.n_rows;
    arma::mat A = S;                              ///< working copy of S
    eigvec.eye(p, p);                             ///< accumulate rotations; starts as identity

    for (int sweep = 0; sweep < 50 * static_cast<int>(p); ++sweep) {
        // Check off-diagonal convergence
        double off = 0.0;
        for (arma::uword q = 0; q < p - 1; ++q)
            for (arma::uword r = q + 1; r < p; ++r)
                off += A(q, r) * A(q, r);
        if (std::sqrt(off) < 1e-13) break;
        detail::jacobi_sweep(A, eigvec);
    }
    eigval = A.diag();

    // Sort ascending
    arma::uvec idx = arma::sort_index(eigval);
    eigval  = eigval(idx);
    eigvec  = eigvec.cols(idx);
#endif
}

/// @}

/// @name QR Factorization
/// @{

// ── qr_full ───────────────────────────────────────────────────────────────────

/// @brief Full (square) QR factorization of an n×p matrix.
///
/// Produces a square n×n orthogonal factor @p Q and an n×p upper-trapezoidal
/// factor @p R such that @c A = Q*R.  This matches the convention of
/// @c arma::qr(), which returns the full (non-economy) @p Q.
///
/// @param[out] Q  n×n orthogonal matrix.
/// @param[out] R  n×p upper-trapezoidal matrix.
/// @param[in]  A  n×p input matrix (not modified).
///
/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::qr directly.
///          Without LAPACK the fallback constructs the full ONB by seeding
///          columns beyond @c A.n_cols with standard basis vectors and then
///          applying modified Gram-Schmidt; numerical quality degrades for
///          nearly rank-deficient inputs.
///
/// @note The caller (@c orthogonal_svd) depends on @p Q being n×n so that
///       @c Q.cols(k, n-1) spans the null space of @c A^T.
inline void qr_full(arma::mat& Q, arma::mat& R, const arma::mat& A) {
#if defined(ARMA_USE_LAPACK)
    arma::qr(Q, R, A);
#else
    // Build a full n×n orthogonal matrix.
    // Step 1: thin QR via modified Gram-Schmidt (gives first k=A.n_cols columns).
    // Step 2: extend to a full ONB using random vectors + re-orthogonalisation.
    const arma::uword n = A.n_rows, p = A.n_cols;

    // Start with all n columns of A padded with random vectors if n > p.
    arma::mat Qfull(n, n, arma::fill::zeros);
    // Copy A's columns first
    for (arma::uword j = 0; j < p && j < n; ++j)
        Qfull.col(j) = A.col(j);
    // Fill remaining columns with standard basis vectors (or random fallback)
    for (arma::uword j = p; j < n; ++j) {
        Qfull(j, j) = 1.0;   ///< standard basis vector e_j used as seed
    }
    // Full MGS orthonormalization
    detail::mgs(Qfull);
    Q = Qfull;

    // R is n×p (upper trapezoidal); only first p rows are non-zero
    R = Q.t() * A;
    // Zero the strictly lower-triangular part of R[0:p, 0:p]
    for (arma::uword j = 0; j < p; ++j)
        for (arma::uword i = j + 1; i < n; ++i)
            R(i, j) = 0.0;
#endif
}

/// @}

/// @name Singular Value Decomposition
/// @{

// ── svd ───────────────────────────────────────────────────────────────────────

/// @brief Full SVD: @c A ≈ U * diag(s) * V', singular values descending.
///
/// Decomposes the n×p matrix @p A into left singular vectors @p U, singular
/// values @p s, and right singular vectors @p V.
///
/// @param[out] U  n×min(n,p) matrix of left singular vectors.
/// @param[out] s  min(n,p)×1 vector of singular values, descending.
/// @param[out] V  p×p matrix of right singular vectors (all p columns).
/// @param[in]  A  n×p input matrix (not modified).
///
/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::svd directly.
///          Without LAPACK the fallback computes @c A^T*A, calls @c eig_sym,
///          and derives @p U via @c U[:,j] = A*V[:,j]/s_j.  Columns with
///          @c s_j < 1e-14 are zeroed.  Accuracy is governed by the
///          conditioning of @c A^T*A (squares the condition number of @p A).
///
/// @note Fallback algorithm: eigendecomposition of \f$A^T A\f$ via cyclic
///       Jacobi (Demmel & Veselić 1992).  Suitable for the small, moderately
///       conditioned matrices typical in ENA (p ≤ ~100).
inline void svd(arma::mat& U, arma::vec& s, arma::mat& V, const arma::mat& A) {
#if defined(ARMA_USE_LAPACK)
    arma::svd(U, s, V, A);
#else
    // One-sided Jacobi on A^T:
    // We work with B = A^T (p×n). The right singular vectors of A are the
    // left singular vectors of B (= columns of U_B after Jacobi on B^T B).
    // We use the eigendecomposition of A^T * A to get V and s.
    const arma::uword n = A.n_rows, p = A.n_cols;
    const arma::uword k = std::min(n, p);

    arma::mat AtA = A.t() * A;   ///< p×p symmetric Gram matrix
    arma::vec eigval;
    arma::mat eigvec;
    eig_sym(eigval, eigvec, AtA);

    // Eigenvalues from eig_sym are ascending; we want descending.
    arma::uvec idx = arma::sort_index(eigval, "descend");
    eigval = eigval(idx);
    eigvec = eigvec.cols(idx);

    // Singular values (non-negative square roots)
    s.set_size(k);
    for (arma::uword j = 0; j < k; ++j)
        s(j) = std::sqrt(std::max(0.0, eigval(j)));

    // V = eigvec (right singular vectors of A are eigvectors of A^T*A)
    V = eigvec;

    // U = A * V / s (left singular vectors, only first k columns needed)
    U.set_size(n, k);
    for (arma::uword j = 0; j < k; ++j) {
        if (s(j) > 1e-14)
            U.col(j) = A * V.col(j) / s(j);
        else
            U.col(j).zeros();
    }
#endif
}

// ── svd_right ─────────────────────────────────────────────────────────────────

/// @brief Singular values and right singular vectors only (left vectors skipped).
///
/// A lighter-weight variant of @c svd() for callers that only need @p s and
/// @p V.  Under LAPACK the full @c arma::svd is still called internally; the
/// savings come in WASM builds where computing @p U is avoided.
///
/// @param[out] s  min(n,p)×1 vector of singular values, descending.
/// @param[out] V  p×p matrix of right singular vectors.
/// @param[in]  A  n×p input matrix (not modified).
///
/// @warning With @c ARMA_USE_LAPACK defined this still invokes the full
///          @c arma::svd; the left-vector output is discarded.
inline void svd_right(arma::vec& s, arma::mat& V, const arma::mat& A) {
#if defined(ARMA_USE_LAPACK)
    arma::mat U;
    arma::svd(U, s, V, A);
#else
    arma::mat U_unused;
    svd(U_unused, s, V, A);
#endif
}

/// @}

/// @name Linear System Solve
/// @{

// ── solve_spd ─────────────────────────────────────────────────────────────────

/// @brief Solve @c A*x = b for symmetric positive (semi-)definite @p A.
///
/// Returns the n×k solution matrix @p x given the n×n system matrix @p A and
/// the n×k right-hand side @p b.  Near-zero eigenvalues are handled gracefully,
/// making this a pseudo-inverse solve safe for rank-deficient Gram matrices.
///
/// @param[in] A    Real symmetric positive (semi-)definite n×n matrix.
/// @param[in] b    n×k right-hand side matrix.
/// @param[in] tol  Relative threshold below which eigenvalues are treated as
///                 zero.  Eigenvalues with @c |λ| < tol * max(|λ|) are
///                 inverted as zero (i.e. those directions are projected out).
///                 Default: 1e-12.
/// @returns  n×k solution matrix @p x satisfying @c A*x ≈ b.
///
/// @warning With @c ARMA_USE_LAPACK defined this delegates to
///          @c arma::solve(A, b, arma::solve_opts::fast), which uses a
///          Cholesky factorization and does **not** handle rank deficiency.
///          Without LAPACK the eigendecomposition-based pseudo-inverse is used,
///          which is safe for rank-deficient inputs but slower.
///
/// @note Without LAPACK: @c A = V * diag(λ) * V^T is computed via @c eig_sym;
///       the solution is @c x = V * diag(1/λ) * V^T * b with small λ skipped.
inline arma::mat solve_spd(const arma::mat& A, const arma::mat& b,
                            double tol = 1e-12) {
#if defined(ARMA_USE_LAPACK)
    return arma::solve(A, b, arma::solve_opts::fast);
#else
    arma::vec  eigval;
    arma::mat  eigvec;
    eig_sym(eigval, eigvec, A);                     ///< A = eigvec * diag(eigval) * eigvec^T

    const double thresh = tol * std::abs(eigval.max());
    arma::vec inv_eigval(eigval.n_elem, arma::fill::zeros);
    for (arma::uword i = 0; i < eigval.n_elem; ++i)
        if (std::abs(eigval(i)) > thresh)
            inv_eigval(i) = 1.0 / eigval(i);

    // x = eigvec * diag(inv_eigval) * eigvec^T * b
    return eigvec * arma::diagmat(inv_eigval) * eigvec.t() * b;
#endif
}

// ── inv_sym ───────────────────────────────────────────────────────────────────

/// @brief Inverse of a symmetric positive (semi-)definite matrix, falling back
///        to the pseudo-inverse when it is singular.
///
/// @param[in] A    Real symmetric positive (semi-)definite n×n matrix.
/// @param[in] tol  Without LAPACK: relative eigenvalue threshold, as in
///                 @c solve_spd.  Default: 1e-12.
/// @returns  n×n matrix: @c A^-1, or the pseudo-inverse when @p A is singular.
///
/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::inv and, when
///          that reports a singular matrix, @c arma::pinv.  Without LAPACK it
///          is the eigendecomposition-based pseudo-inverse (@c solve_spd with
///          the identity), which equals the inverse for a non-singular @p A.
inline arma::mat inv_sym(const arma::mat& A, double tol = 1e-12) {
#if defined(ARMA_USE_LAPACK)
    arma::mat out;
    if (!arma::inv(out, A))
        out = arma::pinv(A);
    return out;
#else
    return solve_spd(A, arma::eye(A.n_rows, A.n_rows), tol);
#endif
}

/// @}

} // namespace linalg
} // namespace qe

#endif // LIBQE_LINALG_FALLBACK_HPP

/**
 * @file normalization.hpp
 * @brief Row-wise normalization utilities for ENA connection networks.
 */
#ifndef LIBQE_NORMALIZATION_HPP
#define LIBQE_NORMALIZATION_HPP

#include <armadillo>
#include <cmath>

namespace qe {

/**
 * @brief Normalize each row to unit L2 norm (project onto the unit hypersphere).
 *
 * Rows whose L2 norm is zero are left unchanged.
 *
 * @param m Input matrix (n_units × n_connections), row-major.
 * @returns Matrix of the same shape with each row divided by its L2 norm.
 * @note Equivalent to @c fun_sphere_norm() in rENA/ena.cpp.
 */
inline arma::mat normalize_networks(arma::mat m) {
    arma::mat out(m.n_rows, m.n_cols, arma::fill::zeros);
    for (arma::uword r = 0; r < m.n_rows; r++) {
        double len = arma::norm(m.row(r), 2);
        if (std::isfinite(len) && len > 0.0) out.row(r) = m.row(r) / len;
    }
    return out;
}

/**
 * @brief Scale every entry by the reciprocal of the largest row L2 norm.
 *
 * Preserves relative magnitudes across rows; does not normalise each row
 * individually.  If the largest row norm is zero, the matrix is returned
 * unchanged.
 *
 * @param m Input matrix (n_units × n_connections), modified in place.
 * @returns Scaled matrix (same shape as @p m).
 * @note Equivalent to @c fun_skip_sphere_norm() in rENA/ena.cpp.
 */
inline arma::mat scale_networks(arma::mat m) {
    double largest = 0.0;
    for (arma::uword r = 0; r < m.n_rows; r++) {
        double len = arma::norm(m.row(r), 2);
        if (std::isfinite(len)) largest = std::max(largest, len);
    }
    if (largest > 0.0) m /= largest;
    return m;
}

} // namespace qe

#endif // LIBQE_NORMALIZATION_HPP

/** @file stability.hpp
 *  @brief Stability metrics and distance-distance matrix correlation for parameter tuning.
 */
#ifndef LIBQE_STABILITY_HPP
#define LIBQE_STABILITY_HPP

#include <armadillo>
#include <cmath>
#include <algorithm>

namespace qe {

/** @brief Compute Pearson correlation between pairwise Euclidean distance structures of two matrices.
 *
 *  Evaluates dist-dist correlation: r(dist(X), dist(Y)), measuring whether relative
 *  neighborhoods and distances are preserved across parameter changes (e.g. window vs lookback).
 *
 *  @param[in] X Numeric matrix (n_rows × n_cols).
 *  @param[in] Y Numeric matrix (n_rows × n_cols).
 *
 *  @returns Pearson correlation coefficient in [-1, 1].
 */
inline double dist_dist_correlation(const arma::mat& X, const arma::mat& Y) {
    size_t n = X.n_rows;
    if (n != Y.n_rows || n < 3) return 0.0;

    size_t n_pairs = n * (n - 1) / 2;
    arma::vec dx(n_pairs);
    arma::vec dy(n_pairs);

    size_t idx = 0;
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = i + 1; j < n; ++j) {
            arma::rowvec diff_x = X.row(i) - X.row(j);
            arma::rowvec diff_y = Y.row(i) - Y.row(j);
            dx[idx] = std::sqrt(arma::dot(diff_x, diff_x));
            dy[idx] = std::sqrt(arma::dot(diff_y, diff_y));
            idx++;
        }
    }

    double mean_x = arma::mean(dx);
    double mean_y = arma::mean(dy);

    arma::vec zx = dx - mean_x;
    arma::vec zy = dy - mean_y;

    double sxx = arma::dot(zx, zx);
    double syy = arma::dot(zy, zy);
    double sxy = arma::dot(zx, zy);

    if (sxx <= 1e-12 || syy <= 1e-12) return 0.0;
    return sxy / std::sqrt(sxx * syy);
}

} // namespace qe

#endif // LIBQE_STABILITY_HPP

/**
 * @file stats.hpp
 * @brief Generic statistics: centering, normal/t quantiles, group confidence
 *        intervals (t-based and IQR outlier fences), and two-group comparison
 *        tests (Welch t-test, Wilcoxon rank-sum).
 *
 * Split out of the former modeling.hpp; the ENA node-position solvers and
 * ena_correlation() are in rENA's libena (positions.hpp).
 */
#ifndef LIBQE_STATS_HPP
#define LIBQE_STATS_HPP

#include <armadillo>
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <vector>

namespace qe {

// ---------------------------------------------------------------------------
/// @name Centering
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Subtract column means so that every dimension is centred at the origin.
 *
 * @param[in] values  Matrix of ENA unit points (n_units × n_dims).
 * @returns A copy of @p values with each column mean subtracted (centre-to-origin).
 *
 * @note Equivalent to `center_data_c()` in rENA/ena.cpp.
 */
inline arma::mat center_points(arma::mat values) {
    return values.each_row() - arma::mean(values);
}

/// @}

/// @cond INTERNAL

// ---------------------------------------------------------------------------
// Normal quantile (probit) — pure C++, no R dependency
// Rational approximation by Peter Acklam; max abs error < 1.15e-9.
// Used by libena's ena_correlation() (rENA) so it doesn't need Rcpp::qnorm.
// ---------------------------------------------------------------------------

inline double normal_quantile(double p) {
    static const double a[] = {
        -3.969683028665376e+01,  2.209460984245205e+02,
        -2.759285104469687e+02,  1.383577518672690e+02,
        -3.066479806614716e+01,  2.506628277459239e+00
    };
    static const double b[] = {
        -5.447609879822406e+01,  1.615858368580409e+02,
        -1.556989798598866e+02,  6.680131188771972e+01,
        -1.328068155288572e+01
    };
    static const double c[] = {
        -7.784894002430293e-03, -3.223964580411365e-01,
        -2.400758277161838e+00, -2.549732539343734e+00,
         4.374664141464968e+00,  2.938163982698783e+00
    };
    static const double d[] = {
         7.784695709041462e-03,  3.224671290700398e-01,
         2.445134137142996e+00,  3.754408661907416e+00
    };
    const double p_low  = 0.02425;
    const double p_high = 1.0 - p_low;
    double q, r;
    if (p < p_low) {
        q = std::sqrt(-2.0 * std::log(p));
        return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
               ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
    } else if (p <= p_high) {
        q = p - 0.5; r = q * q;
        return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
               (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1.0);
    } else {
        q = std::sqrt(-2.0 * std::log(1.0 - p));
        return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
                ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
    }
}

// ---------------------------------------------------------------------------
// Internal helpers: regularized incomplete beta + t-distribution quantile
// ---------------------------------------------------------------------------
//
// These live in qe::detail so they don't pollute the public qe namespace.
// They are used by mean_ci() and are not part of the public API.

namespace detail {

// Lentz continued-fraction evaluation for the regularized incomplete beta.
// Ported from Numerical Recipes in C, §6.4.
inline double betacf(double a, double b, double x) {
    const int    MAXIT = 200;
    const double EPS   = std::numeric_limits<double>::epsilon();
    const double FPMIN = std::numeric_limits<double>::min() / EPS;
    double qab = a + b, qap = a + 1.0, qam = a - 1.0;
    double c = 1.0;
    double d = 1.0 - qab * x / qap;
    if (std::abs(d) < FPMIN) d = FPMIN;
    d = 1.0 / d;
    double h = d;
    for (int m = 1; m <= MAXIT; ++m) {
        const int m2 = 2 * m;
        // even step
        double aa = static_cast<double>(m) * (b - static_cast<double>(m)) * x
                    / ((qam + m2) * (a + m2));
        d = 1.0 + aa * d; if (std::abs(d) < FPMIN) d = FPMIN;
        c = 1.0 + aa / c; if (std::abs(c) < FPMIN) c = FPMIN;
        d = 1.0 / d; h *= d * c;
        // odd step
        aa = -(a + static_cast<double>(m)) * (qab + static_cast<double>(m)) * x
             / ((a + m2) * (qap + m2));
        d = 1.0 + aa * d; if (std::abs(d) < FPMIN) d = FPMIN;
        c = 1.0 + aa / c; if (std::abs(c) < FPMIN) c = FPMIN;
        d = 1.0 / d;
        const double del = d * c;
        h *= del;
        if (std::abs(del - 1.0) <= EPS) break;
    }
    return h;
}

// Regularized incomplete beta I_x(a, b).
inline double betai(double a, double b, double x) {
    if (x <= 0.0) return 0.0;
    if (x >= 1.0) return 1.0;
    const double lbab = std::lgamma(a + b) - std::lgamma(a) - std::lgamma(b);
    const double bt   = std::exp(lbab + a * std::log(x) + b * std::log(1.0 - x));
    if (x < (a + 1.0) / (a + b + 2.0))
        return bt * betacf(a, b, x) / a;
    return 1.0 - bt * betacf(b, a, 1.0 - x) / b;
}

// Quantile function (inverse CDF) for the t-distribution with `df` degrees
// of freedom.  Uses the relation:
//
//   P(T ≤ t | df)  =  1 - I_x(df/2, 1/2) / 2,   x = df / (df + t^2)
//
// Inverted via Newton–Raphson on I_x with bisection fallback; normal-quantile
// plus one Cornish–Fisher correction term provides the starting guess.
inline double t_quantile(double p, double df) {
    const double INF = std::numeric_limits<double>::infinity();
    if (p <= 0.0) return -INF;
    if (p >= 1.0) return  INF;
    if (p == 0.5) return  0.0;
    if (df <= 0.0) return INF;   // df = 0 → t-distribution undefined → treat as ∞

    // Work with pp > 0.5 for numerical stability; negate at the end if needed.
    const bool   flip = (p < 0.5);
    const double pp   = flip ? 1.0 - p : p;
    const double a    = 0.5 * df;
    const double bv   = 0.5;

    // We need x in (0,1) such that I_x(a, bv) = 2*(1-pp).
    const double target = 2.0 * (1.0 - pp);

    // Initial guess via normal quantile + single Cornish–Fisher correction.
    const double z  = normal_quantile(pp);
    const double t0 = z + (z * z * z + z) / (4.0 * df);
    double x = df / (df + t0 * t0);
    x = std::max(1e-12, std::min(1.0 - 1e-12, x));

    // Newton–Raphson with a bisection bracket.
    // f(x) = I_x(a,bv) - target;  f'(x) = x^(a-1)*(1-x)^(bv-1) / B(a,bv)
    const double lbab = std::lgamma(a + bv) - std::lgamma(a) - std::lgamma(bv);
    double xlo = 0.0, xhi = 1.0;
    for (int iter = 0; iter < 100; ++iter) {
        const double fx = betai(a, bv, x) - target;
        // betai is increasing in x, so fx < 0 → x too small → raise lower bound
        if (fx < 0.0) xlo = x; else xhi = x;

        // Derivative: may underflow near boundaries — fall back to bisection.
        const double log_fpx = lbab + (a - 1.0) * std::log(x)
                                     + (bv - 1.0) * std::log(1.0 - x);
        const double fpx = (log_fpx > -700.0) ? std::exp(log_fpx) : 0.0;

        double x_new;
        if (fpx > 0.0) {
            x_new = x - fx / fpx;
            if (x_new <= xlo || x_new >= xhi)
                x_new = 0.5 * (xlo + xhi);
        } else {
            x_new = 0.5 * (xlo + xhi);
        }

        if (std::abs(x_new - x) < 1e-13 * x) break;
        x = x_new;
    }

    const double t_val = std::sqrt(df * (1.0 - x) / x);
    return flip ? -t_val : t_val;
}

// ---------------------------------------------------------------------------
// IQR helper — matches R's quantile(type = 7) (Hyndman–Fan #7), which is
// identical to numpy's np.percentile(method="linear").
//
// For a sorted n-element vector and probability p:
//   h    = (n - 1) * p
//   lo   = floor(h),  frac = h - lo
//   Q(p) = sorted[lo] * (1 - frac) + sorted[lo+1] * frac
// ---------------------------------------------------------------------------

inline double quantile_type7(const arma::vec& sorted, double p) {
    const int n = static_cast<int>(sorted.n_elem);
    if (n == 0) return std::numeric_limits<double>::quiet_NaN();
    if (n == 1) return sorted[0];
    const double h    = (n - 1) * p;
    const int    lo   = static_cast<int>(std::floor(h));
    const double frac = h - lo;
    if (lo + 1 >= n) return sorted[n - 1];          // guard: p == 1.0
    return sorted[lo] * (1.0 - frac) + sorted[lo + 1] * frac;
}

// Interquartile range for a column vector — uses the same type-7 quantile as R.
inline double iqr(const arma::vec& col) {
    const arma::vec s = arma::sort(col);
    return quantile_type7(s, 0.75) - quantile_type7(s, 0.25);
}

} // namespace detail

/// @endcond

// ---------------------------------------------------------------------------
/// @name Group confidence intervals
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Student-t confidence interval for the mean of a group of ENA unit points.
 *
 * For each dimension d:
 * @code
 *   mean ± t_{α/2, n-1} × (sample_sd / sqrt(n))
 * @endcode
 * where α = 1 − @p conf_level and degrees-of-freedom = n − 1.
 *
 * @param[in] points      ENA unit points for a single group, n_units × n_dims.
 * @param[in] conf_level  Confidence level for the interval (default 0.95).
 * @returns An n_dims × 3 matrix whose columns are [mean, ci_lower, ci_upper].
 *          When n == 0 all entries are NaN.
 *          When n == 1 the CI bounds are ±Inf.
 */
inline arma::mat mean_ci(const arma::mat& points, double conf_level = 0.95) {
    const int n      = static_cast<int>(points.n_rows);
    const int n_dims = static_cast<int>(points.n_cols);

    arma::mat out(n_dims, 3);
    out.fill(std::numeric_limits<double>::quiet_NaN());

    if (n == 0) return out;

    const double df     = static_cast<double>(n - 1);
    const double t_crit = detail::t_quantile(1.0 - (1.0 - conf_level) / 2.0, df);
    const double INF    = std::numeric_limits<double>::infinity();

    for (int d = 0; d < n_dims; ++d) {
        const arma::vec col = points.col(d);
        const double mu = arma::mean(col);
        out(d, 0) = mu;
        if (n == 1) {
            // t_crit = Inf and se = 0 would give NaN; explicitly set ±Inf
            out(d, 1) = -INF;
            out(d, 2) =  INF;
        } else {
            const double se = arma::stddev(col) / std::sqrt(static_cast<double>(n));
            out(d, 1) = mu - t_crit * se;
            out(d, 2) = mu + t_crit * se;
        }
    }
    return out;
}

/**
 * @brief Tukey-fence outlier interval for a group of ENA unit points.
 *
 * For each dimension d the half-width is computed as:
 * @code
 *   half_width[d] = iqr_factor * IQR(points[:, d])
 * @endcode
 * and the interval [−half_width, +half_width] is symmetric around zero,
 * matching rENA's formula:
 * @code
 *   outlier.interval.values = matrix(...) * c(-1, 1)
 * @endcode
 * where the matrix is built from `c(IQR(dim1), IQR(dim2), ...) * iqr_factor`.
 *
 * IQR uses R's default type-7 / Hyndman-Fan #7 quantile, identical to
 * `numpy.percentile(method="linear")`.
 *
 * @param[in] points      ENA unit points for a single group, n_units × n_dims.
 * @param[in] iqr_factor  Multiplier applied to the IQR (default 1.5, the
 *                        standard Tukey fence threshold).
 * @returns An n_dims × 2 matrix whose columns are [lower, upper].
 *          When n == 0 all entries are NaN.
 *          When n == 1, IQR == 0, so lower == upper == 0.
 */
inline arma::mat outlier_ci(const arma::mat& points, double iqr_factor = 1.5) {
    const int n_dims = static_cast<int>(points.n_cols);
    arma::mat out(n_dims, 2, arma::fill::zeros);

    if (points.n_rows == 0) {
        out.fill(std::numeric_limits<double>::quiet_NaN());
        return out;
    }

    for (int d = 0; d < n_dims; ++d) {
        const double half_width = iqr_factor * detail::iqr(points.col(d));
        out(d, 0) = -half_width;
        out(d, 1) =  half_width;
    }
    return out;
}

/// @}

// ---------------------------------------------------------------------------
/// @name Two-group comparison statistics
/// @{
// ---------------------------------------------------------------------------

/**
 * @brief Per-dimension statistics from a two-group comparison.
 *
 * Matches the output of rENA-api's `group.stats()`:
 *
 * - **Parametric** (Welch two-sample t-test): `t`, `df`, `pvalue_t`,
 *   `cohens_d` (pooled-SD formula), `means`, `sds`.
 * - **Non-parametric** (Wilcoxon rank-sum / Mann-Whitney U):
 *   `U` (= R's W statistic for group 1), `pvalue_u` (normal approximation
 *   with tie correction and continuity correction), `effect_r`
 *   (rank-biserial: `1 − 2*U / (n1*n2)`), `medians`.
 *
 * All per-dimension vectors have length `n_dims`; `means`, `sds`, and
 * `medians` are `2 × n_dims` matrices with row 0 = group 1, row 1 = group 2.
 * Entries are NaN when a statistic is undefined (e.g. n < 2 for parametric
 * tests, or n = 0 for either group).
 */
struct GroupStats {
    int n1;       ///< Sample size of group 1.
    int n2;       ///< Sample size of group 2.

    // Parametric
    arma::vec t;          ///< Welch two-sample t-statistics.
    arma::vec df;         ///< Welch–Satterthwaite degrees of freedom.
    arma::vec pvalue_t;   ///< Two-tailed p-values from the Welch t-test.
    arma::vec cohens_d;   ///< Cohen's d (pooled-SD formula, group1 − group2).
    arma::mat means;      ///< 2 × n_dims: means; row 0 = group1, row 1 = group2.
    arma::mat sds;        ///< 2 × n_dims: sample standard deviations.

    // Non-parametric
    arma::vec U;          ///< Wilcoxon rank-sum U for group 1 (= R's W statistic).
    arma::vec pvalue_u;   ///< Two-tailed p-values (normal approx, tie + continuity correction).
    arma::vec effect_r;   ///< Rank-biserial effect: 1 − 2·U / (n1·n2).
    arma::mat medians;    ///< 2 × n_dims: medians; row 0 = group1, row 1 = group2.
};

/// @cond INTERNAL
namespace detail {

// Two-tailed p-value for a Welch t-statistic.
// Uses the relation P(T≤t|df) = 1 − I_x(df/2, 1/2)/2, x = df/(df+t²).
inline double t_test_pvalue(double t_stat, double df) {
    if (!std::isfinite(t_stat) || df <= 0.0)
        return std::numeric_limits<double>::quiet_NaN();
    const double x = df / (df + t_stat * t_stat);
    return betai(df / 2.0, 0.5, x);   // == 2 * P(T_df > |t_stat|)
}

// Wilcoxon rank-sum U (= R's W for group 1) and two-tailed p-value.
// Normal approximation with tie correction and continuity correction,
// matching R's wilcox.test(..., exact=FALSE, correct=TRUE).
struct WilcoxResult { double U; double p; };

inline WilcoxResult wilcox_ranksum(const arma::vec& g1, const arma::vec& g2) {
    const int    n1  = static_cast<int>(g1.n_elem);
    const int    n2  = static_cast<int>(g2.n_elem);
    const int    n   = n1 + n2;
    const double NaN = std::numeric_limits<double>::quiet_NaN();

    if (n1 == 0 || n2 == 0) return {NaN, NaN};

    // Merge both groups into a single array, tagging each element's origin.
    arma::vec      vals(n);
    std::vector<int> grp(n);
    for (int i = 0; i < n1; ++i) { vals[i]      = g1[i]; grp[i]      = 0; }
    for (int i = 0; i < n2; ++i) { vals[n1 + i] = g2[i]; grp[n1 + i] = 1; }

    arma::uvec ord = arma::sort_index(vals);

    // Assign average ranks for tied values; accumulate tie correction T = Σ(t³−t).
    arma::vec rnk(n);
    double    T_ties = 0.0;
    arma::uword i = 0;
    while (i < static_cast<arma::uword>(n)) {
        arma::uword j = i;
        while (j + 1 < static_cast<arma::uword>(n) &&
               vals[ord[j + 1]] == vals[ord[j]])
            ++j;
        const double avg = 0.5 * (static_cast<double>(i) + static_cast<double>(j)) + 1.0;
        for (arma::uword k = i; k <= j; ++k) rnk[ord[k]] = avg;
        const double t = static_cast<double>(j - i + 1);
        T_ties += t * (t * t - 1.0);   // t³ − t
        i = j + 1;
    }

    // W = rank-sum for group 1; U = W − n1*(n1+1)/2 (R's wilcox.test statistic).
    double W = 0.0;
    for (int k = 0; k < n; ++k)
        if (grp[k] == 0) W += rnk[k];
    const double U = W - static_cast<double>(n1) * (n1 + 1) / 2.0;

    // Normal approximation: Var[U] = n1·n2/12 · ((n+1) − T_ties/(n·(n−1)))
    const double mu_U  = static_cast<double>(n1) * n2 / 2.0;
    const double tie_adj = (n > 1) ? T_ties / (static_cast<double>(n) * (n - 1)) : 0.0;
    double var_U = static_cast<double>(n1) * n2 / 12.0
                   * (static_cast<double>(n + 1) - tie_adj);
    if (var_U < 0.0) var_U = 0.0;

    double p;
    if (var_U == 0.0) {
        p = (U == mu_U) ? 1.0 : 0.0;
    } else {
        // Continuity correction: |U − mu| − 0.5
        double z = (std::abs(U - mu_U) - 0.5) / std::sqrt(var_U);
        if (z < 0.0) z = 0.0;
        p = std::erfc(z / std::sqrt(2.0));   // = 2*(1−Φ(z))
        p = std::min(p, 1.0);
    }

    return {U, p};
}

} // namespace detail (group stats helpers)
/// @endcond

/**
 * @brief Compute per-dimension two-group comparison statistics.
 *
 * For each dimension @p d of the input matrices the function computes:
 *
 * **Parametric (Welch two-sample t-test)**
 *  - t-statistic and Welch–Satterthwaite degrees of freedom
 *  - Two-tailed p-value via the regularised incomplete beta function
 *  - Cohen's d using the pooled-SD formula
 *    (`(μ1−μ2) / √(((n1−1)s1² + (n2−1)s2²)/(n1+n2−2))`)
 *  - Per-group means and sample standard deviations
 *
 * **Non-parametric (Wilcoxon rank-sum / Mann-Whitney U)**
 *  - U statistic for group 1 (identical to R's `wilcox.test` W statistic)
 *  - Two-tailed p-value via the normal approximation with tie correction
 *    and continuity correction (matching R's `wilcox.test(..., correct=TRUE)`)
 *  - Rank-biserial effect size: `1 − 2·U / (n1·n2)`
 *    (rENA-api's `nonparam.effect` formula)
 *  - Per-group medians
 *
 * @param[in] g1  Unit points for group 1, n1 × n_dims.
 * @param[in] g2  Unit points for group 2, n2 × n_dims.
 * @returns A @ref GroupStats struct.
 *
 * @note Parametric statistics are NaN when either group has fewer than 2 units
 *       (no sample variance).  All statistics are NaN when either group is empty.
 * @note The Wilcoxon p-value uses the normal approximation regardless of sample
 *       size; R switches to an exact test for small (n ≤ ~50) tie-free samples,
 *       so p-values may differ slightly in those cases.
 * @note Equivalent to `group.stats()` in rENA-api, moved to the C++ layer so
 *       all language bindings share a single implementation.
 */
inline GroupStats group_stats(const arma::mat& g1, const arma::mat& g2) {
    const int    n1     = static_cast<int>(g1.n_rows);
    const int    n2     = static_cast<int>(g2.n_rows);
    const int    n_dims = static_cast<int>(g1.n_cols);
    const double NaN    = std::numeric_limits<double>::quiet_NaN();

    GroupStats out;
    out.n1 = n1;
    out.n2 = n2;

    out.t        = arma::vec(n_dims); out.t.fill(NaN);
    out.df       = arma::vec(n_dims); out.df.fill(NaN);
    out.pvalue_t = arma::vec(n_dims); out.pvalue_t.fill(NaN);
    out.cohens_d = arma::vec(n_dims); out.cohens_d.fill(NaN);
    out.means    = arma::mat(2, n_dims); out.means.fill(NaN);
    out.sds      = arma::mat(2, n_dims); out.sds.fill(NaN);
    out.U        = arma::vec(n_dims); out.U.fill(NaN);
    out.pvalue_u = arma::vec(n_dims); out.pvalue_u.fill(NaN);
    out.effect_r = arma::vec(n_dims); out.effect_r.fill(NaN);
    out.medians  = arma::mat(2, n_dims); out.medians.fill(NaN);

    if (n1 == 0 || n2 == 0) return out;

    for (int d = 0; d < n_dims; ++d) {
        const arma::vec c1 = g1.col(d);
        const arma::vec c2 = g2.col(d);

        // ── Means and medians ───────────────────────────────────────────────
        const double m1 = arma::mean(c1);
        const double m2 = arma::mean(c2);
        out.means(0, d) = m1;
        out.means(1, d) = m2;
        out.medians(0, d) = arma::median(c1);
        out.medians(1, d) = arma::median(c2);

        // ── Standard deviations (ddof = 1, matching R's sd()) ───────────────
        const double s1 = (n1 > 1) ? arma::stddev(c1) : NaN;
        const double s2 = (n2 > 1) ? arma::stddev(c2) : NaN;
        out.sds(0, d) = s1;
        out.sds(1, d) = s2;

        // ── Welch two-sample t-test ──────────────────────────────────────────
        if (n1 > 1 && n2 > 1 && std::isfinite(s1) && std::isfinite(s2)) {
            const double v1n = s1 * s1 / n1;
            const double v2n = s2 * s2 / n2;
            const double se  = std::sqrt(v1n + v2n);
            if (se > 0.0) {
                const double t_stat = (m1 - m2) / se;
                const double df_val = (v1n + v2n) * (v1n + v2n)
                                    / (v1n * v1n / (n1 - 1) + v2n * v2n / (n2 - 1));
                out.t(d)        = t_stat;
                out.df(d)       = df_val;
                out.pvalue_t(d) = detail::t_test_pvalue(t_stat, df_val);
            }
        }

        // ── Cohen's d (pooled SD) ────────────────────────────────────────────
        if (n1 + n2 > 2 && std::isfinite(s1) && std::isfinite(s2)) {
            const double pooled_var = ((n1 - 1) * s1 * s1 + (n2 - 1) * s2 * s2)
                                    / (n1 + n2 - 2);
            const double pooled_sd  = std::sqrt(pooled_var);
            out.cohens_d(d) = (pooled_sd > 0.0) ? (m1 - m2) / pooled_sd : 0.0;
        }

        // ── Wilcoxon rank-sum ────────────────────────────────────────────────
        const detail::WilcoxResult wr = detail::wilcox_ranksum(c1, c2);
        out.U(d)        = wr.U;
        out.pvalue_u(d) = wr.p;
        if (std::isfinite(wr.U))
            out.effect_r(d) = 1.0 - 2.0 * wr.U / (static_cast<double>(n1) * n2);
    }

    return out;
}

/// @}

} // namespace qe

#endif // LIBQE_STATS_HPP

/** @file trajectory_distance.hpp
 *  @brief Integrated Euclidean distance between trajectory curves and lagged distance metrics.
 */
#ifndef LIBQE_TRAJECTORY_DISTANCE_HPP
#define LIBQE_TRAJECTORY_DISTANCE_HPP

#include <armadillo>
#include <vector>
#include <cmath>
#include <algorithm>
#include "trajectory.hpp"
#include "linalg_fallback.hpp"

namespace qe {

struct GaussLegendreRule {
    arma::vec nodes;
    arma::vec weights;
};

/** @brief Obtain exact 32-point Gauss-Legendre quadrature nodes and weights on [-1, 1] via Golub-Welsch. */
inline const GaussLegendreRule& get_gl_rule_32() {
    static const GaussLegendreRule rule = []() {
        int N = 32;
        arma::mat J(N, N, arma::fill::zeros);
        for (int k = 1; k < N; ++k) {
            double b = static_cast<double>(k) / std::sqrt(4.0 * k * k - 1.0);
            J(k, k - 1) = b;
            J(k - 1, k) = b;
        }
        arma::vec evals;
        arma::mat evecs;
        qe::linalg::eig_sym(evals, evecs, J);   // LAPACK-free in WASM builds
        arma::vec w = 2.0 * arma::square(evecs.row(0).t());
        return GaussLegendreRule{evals, w};
    }();
    return rule;
}

/** @brief Compute integrated Euclidean distance between two 2D parametric polynomial curves.
 *
 *  Computes integral_{t_start}^{t_end} ||F_a(t) - F_b(t)|| dt
 *  using 32-point Gauss-Legendre quadrature.
 *
 *  @param[in] coeffs_ax  Coefficients of x_a(t).
 *  @param[in] coeffs_ay  Coefficients of y_a(t).
 *  @param[in] coeffs_bx  Coefficients of x_b(t).
 *  @param[in] coeffs_by  Coefficients of y_b(t).
 *  @param[in] t_start    Start of integration interval (default 0.0).
 *  @param[in] t_end      End of integration interval (default 1.0).
 *
 *  @returns Integrated distance scalar.
 */
inline double integrated_curve_distance(
    const arma::vec& coeffs_ax,
    const arma::vec& coeffs_ay,
    const arma::vec& coeffs_bx,
    const arma::vec& coeffs_by,
    double t_start = 0.0,
    double t_end = 1.0
) {
    if (t_end <= t_start) return 0.0;

    const auto& rule = get_gl_rule_32();
    double mid = 0.5 * (t_start + t_end);
    double half_len = 0.5 * (t_end - t_start);
    double total = 0.0;

    for (size_t i = 0; i < rule.nodes.n_elem; ++i) {
        double t_eval = mid + half_len * rule.nodes[i];
        double xa = eval_poly_scalar(coeffs_ax, t_eval);
        double ya = eval_poly_scalar(coeffs_ay, t_eval);
        double xb = eval_poly_scalar(coeffs_bx, t_eval);
        double yb = eval_poly_scalar(coeffs_by, t_eval);
        double dx = xa - xb;
        double dy = ya - yb;
        double dist = std::sqrt(dx * dx + dy * dy);

        total += rule.weights[i] * dist;
    }

    return half_len * total;
}

/** @brief Compute lagged curve distance: 1/(1-lag) * integral_{lag}^1 ||F_follower(t) - F_leader(t - lag)|| dt. */
inline double lagged_curve_distance(
    const arma::vec& coeffs_fol_x,
    const arma::vec& coeffs_fol_y,
    const arma::vec& coeffs_ldr_x,
    const arma::vec& coeffs_ldr_y,
    double lag = 0.0
) {
    if (lag < 0.0 || lag >= 0.999) return 0.0;
    double t_start = lag;
    double t_end = 1.0;
    const auto& rule = get_gl_rule_32();
    double mid = 0.5 * (t_start + t_end);
    double half_len = 0.5 * (t_end - t_start);
    double total = 0.0;

    for (size_t i = 0; i < rule.nodes.n_elem; ++i) {
        double t_fol = mid + half_len * rule.nodes[i];
        double t_ldr = t_fol - lag;

        double x_fol = eval_poly_scalar(coeffs_fol_x, t_fol);
        double y_fol = eval_poly_scalar(coeffs_fol_y, t_fol);
        double x_ldr = eval_poly_scalar(coeffs_ldr_x, t_ldr);
        double y_ldr = eval_poly_scalar(coeffs_ldr_y, t_ldr);

        double dx = x_fol - x_ldr;
        double dy = y_fol - y_ldr;
        double dist = std::sqrt(dx * dx + dy * dy);

        total += rule.weights[i] * dist;
    }

    double integral = half_len * total;
    double interval = 1.0 - lag;
    return (interval > 1e-6) ? (integral / interval) : 0.0;
}

/** @brief Compute pairwise integrated distance matrix across K polynomial trajectory curves. */
inline arma::mat pairwise_trajectory_distance_matrix(
    const std::vector<arma::vec>& all_coeffs_x,
    const std::vector<arma::vec>& all_coeffs_y
) {
    size_t K = all_coeffs_x.size();
    arma::mat D(K, K, arma::fill::zeros);

    for (size_t i = 0; i < K; ++i) {
        for (size_t j = i + 1; j < K; ++j) {
            double dist = integrated_curve_distance(
                all_coeffs_x[i], all_coeffs_y[i],
                all_coeffs_x[j], all_coeffs_y[j]
            );
            D(i, j) = dist;
            D(j, i) = dist;
        }
    }

    return D;
}

} // namespace qe

#endif // LIBQE_TRAJECTORY_DISTANCE_HPP

/** @file trajectory_following.hpp
 *  @brief Who-follows-whom leader-follower turn lag analysis along shared timeline.
 */
#ifndef LIBQE_TRAJECTORY_FOLLOWING_HPP
#define LIBQE_TRAJECTORY_FOLLOWING_HPP

#include <armadillo>
#include <vector>
#include <cmath>
#include <algorithm>
#include <limits>

namespace qe {

/** @brief Result of signed turn-lag sweep across delta in [-max_lag, max_lag]. */
struct SignedTurnLagResult {
    int best_lag = 0;             ///< delta* that minimizes mean distance
    double min_mean_distance = 0.0;
    arma::ivec lags;              ///< Range of evaluated lags
    arma::vec mean_distances;     ///< Mean distance for each lag
    arma::ivec valid_counts;      ///< Number of paired turns for each lag
};

/** @brief Compute mean distance between agent A at turn T and agent B at turn T - delta.
 *
 *  @param[in] pts_a    n_a × 2 coordinates for agent A.
 *  @param[in] pts_b    n_b × 2 coordinates for agent B.
 *  @param[in] times_a  Length n_a vector of calendar/turn indices for A.
 *  @param[in] times_b  Length n_b vector of calendar/turn indices for B.
 *  @param[in] delta    Signed lag in turns (delta > 0: A compares to B's past => A follows B;
 *                      delta < 0: B compares to A's past => B follows A; delta = 0: sync).
 *
 *  @returns Pair of (mean_distance, valid_pairs_count).
 */
inline std::pair<double, int> signed_turn_lag_distance(
    const arma::mat& pts_a,
    const arma::mat& pts_b,
    const arma::vec& times_a,
    const arma::vec& times_b,
    int delta
) {
    int n_a = static_cast<int>(pts_a.n_rows);
    int n_b = static_cast<int>(pts_b.n_rows);
    if (n_a == 0 || n_b == 0) return {0.0, 0};

    double sum_dist = 0.0;
    int count = 0;

    for (int i = 0; i < n_a; ++i) {
        double t_a = times_a[i];
        double target_tb = t_a - static_cast<double>(delta);

        // Find the latest turn for agent B that occurred at or before target_tb
        int best_b_idx = -1;
        for (int j = 0; j < n_b; ++j) {
            if (times_b[j] <= target_tb) {
                best_b_idx = j;
            } else {
                break;
            }
        }

        if (best_b_idx >= 0) {
            double dx = pts_a(i, 0) - pts_b(best_b_idx, 0);
            double dy = pts_a(i, 1) - pts_b(best_b_idx, 1);
            sum_dist += std::sqrt(dx * dx + dy * dy);
            count++;
        }
    }

    double mean_d = (count > 0) ? (sum_dist / count) : 0.0;
    return {mean_d, count};
}

/** @brief Sweep signed turn lags in [-max_lag, max_lag] to determine leader-follower dynamics. */
inline SignedTurnLagResult best_signed_turn_lag(
    const arma::mat& pts_a,
    const arma::mat& pts_b,
    const arma::vec& times_a,
    const arma::vec& times_b,
    int max_lag = 15
) {
    int k_max = std::max(1, max_lag);
    int n_lags = 2 * k_max + 1;

    SignedTurnLagResult res;
    res.lags.set_size(n_lags);
    res.mean_distances.set_size(n_lags);
    res.valid_counts.set_size(n_lags);

    double min_dist = std::numeric_limits<double>::infinity();
    int best_delta = 0;

    for (int idx = 0; idx < n_lags; ++idx) {
        int delta = -k_max + idx;
        res.lags[idx] = delta;

        auto [mean_d, count] = signed_turn_lag_distance(pts_a, pts_b, times_a, times_b, delta);
        res.mean_distances[idx] = mean_d;
        res.valid_counts[idx] = count;

        if (count > 0 && mean_d < min_dist) {
            min_dist = mean_d;
            best_delta = delta;
        }
    }

    res.best_lag = best_delta;
    res.min_mean_distance = (min_dist == std::numeric_limits<double>::infinity()) ? 0.0 : min_dist;
    return res;
}

} // namespace qe

#endif // LIBQE_TRAJECTORY_FOLLOWING_HPP

/** @file trajectory.hpp
 *  @brief Parametric trajectory curve fitting, LOOCV-2D degree selection,
 *         and differential geometry (velocity, acceleration, curvature, turns).
 */
#ifndef LIBQE_TRAJECTORY_HPP
#define LIBQE_TRAJECTORY_HPP

#include <armadillo>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>
#include <limits>

#include "linalg_fallback.hpp"

namespace qe {

/** @brief Polynomial curve fit result in 2D epistemic space. */
struct PolyCurveFit {
    int degree = 1;
    arma::vec coeffs_x;     ///< Raw coefficients [c0, ..., cd] for x(t) = sum(ck * t^k)
    arma::vec coeffs_y;     ///< Raw coefficients [c0, ..., cd] for y(t) = sum(ck * t^k)
    std::string basis = "orthogonal"; ///< Basis used during fitting: "orthogonal" (R stats::poly-like) or "raw".
    arma::vec basis_alpha;  ///< Orthogonal-basis recurrence alpha values (empty for raw basis).
    arma::vec basis_norm2;  ///< Orthogonal-basis recurrence squared norms, matching stats::poly coefs$norm2.
    double cv_error = 0.0;  ///< 2D LOOCV error
    double aic = 0.0;       ///< Joint bivariate AIC
    arma::vec t;            ///< Evaluated or input time points
    arma::mat fitted;       ///< n × 2 fitted coordinates
};

inline arma::mat build_raw_poly_design(const arma::vec& t_vec, int degree) {
    int n = static_cast<int>(t_vec.n_elem);
    arma::mat X(n, degree + 1);
    for (int i = 0; i < n; ++i) {
        double cur_t = t_vec[i];
        double p = 1.0;
        for (int k = 0; k <= degree; ++k) {
            X(i, k) = p;
            p *= cur_t;
        }
    }
    return X;
}

inline arma::vec poly_multiply_by_linear(const arma::vec& coeffs, double shift) {
    arma::vec out(coeffs.n_elem + 1, arma::fill::zeros);
    for (size_t k = 0; k < coeffs.n_elem; ++k) {
        out[k] += -shift * coeffs[k];
        out[k + 1] += coeffs[k];
    }
    return out;
}

inline arma::mat orthogonal_poly_basis_from_coefs(
    const arma::vec& t_vec,
    int degree,
    const arma::vec& alpha,
    const arma::vec& norm2
) {
    int n = static_cast<int>(t_vec.n_elem);
    arma::mat Z(n, degree + 1, arma::fill::ones);
    if (degree >= 1) {
        Z.col(1) = t_vec - alpha[0];
    }
    for (int i = 2; i <= degree; ++i) {
        Z.col(i) = (t_vec - alpha[i - 1]) % Z.col(i - 1) -
                   (norm2[i] / norm2[i - 1]) * Z.col(i - 2);
    }
    for (int i = 0; i <= degree; ++i) {
        Z.col(i) /= std::sqrt(norm2[i + 1]);
    }
    return Z.cols(1, degree);
}

inline arma::mat orthogonal_basis_raw_coefficients(
    int degree,
    const arma::vec& alpha,
    const arma::vec& norm2
) {
    arma::mat raw(degree + 1, degree + 1, arma::fill::zeros);
    raw(0, 0) = 1.0;
    if (degree >= 1) {
        arma::vec p1(2);
        p1[0] = -alpha[0];
        p1[1] = 1.0;
        raw.submat(0, 1, 1, 1) = p1 / std::sqrt(norm2[2]);
    }
    arma::vec prev2(1);
    prev2[0] = 1.0;
    arma::vec prev1;
    if (degree >= 1) {
        prev1.set_size(2);
        prev1[0] = -alpha[0];
        prev1[1] = 1.0;
    }
    for (int i = 2; i <= degree; ++i) {
        arma::vec term = poly_multiply_by_linear(prev1, alpha[i - 1]);
        arma::vec padded_prev2(term.n_elem, arma::fill::zeros);
        padded_prev2.subvec(0, prev2.n_elem - 1) = prev2;
        arma::vec current = term - (norm2[i] / norm2[i - 1]) * padded_prev2;
        raw.submat(0, i, i, i) = current / std::sqrt(norm2[i + 1]);
        prev2 = prev1;
        prev1 = current;
    }
    return raw;
}

inline bool modified_gram_schmidt_qr(
    const arma::mat& X,
    arma::mat& Q,
    arma::mat& R
) {
    const int n = static_cast<int>(X.n_rows);
    const int p = static_cast<int>(X.n_cols);
    Q.zeros(n, p);
    R.zeros(p, p);

    for (int k = 0; k < p; ++k) {
        arma::vec v = X.col(k);
        for (int j = 0; j < k; ++j) {
            R(j, k) = arma::dot(Q.col(j), v);
            v -= R(j, k) * Q.col(j);
        }

        double norm = std::sqrt(arma::dot(v, v));
        if (norm <= 1e-12) {
            return false;
        }

        R(k, k) = norm;
        Q.col(k) = v / norm;
    }

    return true;
}

inline arma::mat build_orthogonal_poly_design(
    const arma::vec& t_vec,
    int degree,
    arma::vec& alpha,
    arma::vec& norm2
) {
    int n = static_cast<int>(t_vec.n_elem);
    if (degree < 1) {
        return arma::mat(n, 0);
    }
    arma::vec unique_t = arma::unique(t_vec);
    if (static_cast<size_t>(degree) >= unique_t.n_elem) {
        throw std::invalid_argument("orthogonal polynomial degree must be less than number of unique time points.");
    }

    double xbar = arma::mean(t_vec);
    arma::vec centered = t_vec - xbar;
    arma::mat X = build_raw_poly_design(centered, degree);

    arma::mat Q;
    arma::mat R;
    bool ok = modified_gram_schmidt_qr(X, Q, R);
    if (!ok) {
        throw std::invalid_argument("orthogonal polynomial degree must be less than number of unique time points.");
    }

    arma::mat Z = Q * arma::diagmat(R.diag());
    arma::vec z_norm2 = arma::sum(arma::square(Z), 0).t();

    alpha.set_size(degree);
    for (int i = 0; i < degree; ++i) {
        alpha[i] = arma::dot(centered, arma::square(Z.col(i))) / z_norm2[i] + xbar;
    }

    norm2.set_size(degree + 2);
    norm2[0] = 1.0;
    for (int i = 0; i <= degree; ++i) {
        norm2[i + 1] = z_norm2[i];
    }

    return orthogonal_poly_basis_from_coefs(t_vec, degree, alpha, norm2);
}

inline arma::mat build_poly_design(
    const arma::vec& t_vec,
    int degree,
    const std::string& basis,
    arma::vec& alpha,
    arma::vec& norm2,
    arma::mat& basis_raw_coeffs
) {
    int n = static_cast<int>(t_vec.n_elem);
    arma::mat X(n, degree + 1);
    X.col(0).ones();
    if (basis == "orthogonal") {
        arma::mat Z = build_orthogonal_poly_design(t_vec, degree, alpha, norm2);
        X.cols(1, degree) = Z;
        basis_raw_coeffs = orthogonal_basis_raw_coefficients(degree, alpha, norm2);
    } else if (basis == "raw") {
        X = build_raw_poly_design(t_vec, degree);
        alpha.reset();
        norm2.reset();
        basis_raw_coeffs.eye(degree + 1, degree + 1);
    } else {
        throw std::invalid_argument("basis must be 'orthogonal' or 'raw'.");
    }
    return X;
}

/** @brief Differential geometry profile along a trajectory curve. */
struct TrajectoryDerivatives {
    arma::vec t;
    arma::vec vx;           ///< dx/dt
    arma::vec vy;           ///< dy/dt
    arma::vec speed;        ///< sqrt(vx^2 + vy^2)
    arma::vec ax;           ///< d^2x/dt^2
    arma::vec ay;           ///< d^2y/dt^2
    arma::vec heading_rate; ///< omega = vx*ay - vy*ax (angular turning velocity)
    arma::vec curvature;    ///< kappa = omega / max(speed^3, eps)
};

/** @brief Evaluate a 1D polynomial sum(c_k * t^k) at a single time t. */
inline double eval_poly_scalar(const arma::vec& coeffs, double t) {
    double val = 0.0;
    double t_pow = 1.0;
    for (size_t k = 0; k < coeffs.n_elem; ++k) {
        val += coeffs[k] * t_pow;
        t_pow *= t;
    }
    return val;
}

/** @brief Evaluate polynomial 1st derivative sum(k * c_k * t^(k-1)) at a single time t. */
inline double eval_poly_deriv1_scalar(const arma::vec& coeffs, double t) {
    double val = 0.0;
    double t_pow = 1.0;
    for (size_t k = 1; k < coeffs.n_elem; ++k) {
        val += static_cast<double>(k) * coeffs[k] * t_pow;
        t_pow *= t;
    }
    return val;
}

/** @brief Evaluate polynomial 2nd derivative sum(k*(k-1) * c_k * t^(k-2)) at a single time t. */
inline double eval_poly_deriv2_scalar(const arma::vec& coeffs, double t) {
    double val = 0.0;
    double t_pow = 1.0;
    for (size_t k = 2; k < coeffs.n_elem; ++k) {
        val += static_cast<double>(k * (k - 1)) * coeffs[k] * t_pow;
        t_pow *= t;
    }
    return val;
}

/** @brief Evaluate 2D parametric curve (x(t), y(t)) at multiple time points. */
inline arma::mat eval_poly_curve(
    const arma::vec& coeffs_x,
    const arma::vec& coeffs_y,
    const arma::vec& t_eval
) {
    size_t n = t_eval.n_elem;
    arma::mat out(n, 2);
    for (size_t i = 0; i < n; ++i) {
        double t = t_eval[i];
        out(i, 0) = eval_poly_scalar(coeffs_x, t);
        out(i, 1) = eval_poly_scalar(coeffs_y, t);
    }
    return out;
}

/** @brief Fit a 2D polynomial curve with automatic LOOCV-2D or AIC degree selection.
 *
 *  Given coordinates (x_i, y_i) and times t_i in [0, 1], constructs polynomial models
 *  for degrees d in [1, max_degree]. LOOCV error is computed via the closed-form hat-matrix formula:
 *  CV(d) = sum( (r_{x,i} / (1 - h_{ii}))^2 + (r_{y,i} / (1 - h_{ii}))^2 )
 *
 *  @param[in] points      n × 2 matrix of trajectory coordinates.
 *  @param[in] t           Vector of length n of strictly increasing times (default [0, 1] linearly spaced).
 *  @param[in] max_degree  Maximum polynomial degree to evaluate (default 3).
 *  @param[in] fixed_degree When >= 1, uses this exact degree without search.
 *  @param[in] criterion   "loocv" (default) or "aic".
 *
 *  @returns PolyCurveFit containing optimal degree, coefficients, errors, and fitted points.
 */
inline PolyCurveFit fit_poly_loocv2d(
    const arma::mat& points,
    const arma::vec& t = arma::vec(),
    int max_degree = 3,
    int fixed_degree = 0,
    const std::string& criterion = "loocv",
    const std::string& basis = "orthogonal"
) {
    int n = static_cast<int>(points.n_rows);
    if (n < 2 || points.n_cols < 2) {
        throw std::invalid_argument("fit_poly_loocv2d requires at least 2 points and 2 columns.");
    }

    // Default normalized time t in [0, 1]
    arma::vec t_vec;
    if (t.n_elem == static_cast<size_t>(n)) {
        t_vec = t;
    } else {
        t_vec = arma::linspace<arma::vec>(0.0, 1.0, n);
    }

    arma::vec x = points.col(0);
    arma::vec y = points.col(1);

    int deg_upper = std::min(max_degree, std::max(1, n - 2));
    if (fixed_degree >= 1) {
        deg_upper = std::min(fixed_degree, n - 1);
    }

    int best_deg = 1;
    double best_score = std::numeric_limits<double>::infinity();
    PolyCurveFit best_fit;

    int deg_start = (fixed_degree >= 1) ? fixed_degree : 1;

    for (int d = deg_start; d <= deg_upper; ++d) {
        int n_params = d + 1;
        if (n_params > n) break;

        arma::vec basis_alpha;
        arma::vec basis_norm2;
        arma::mat basis_raw_coeffs;
        arma::mat X = build_poly_design(t_vec, d, basis, basis_alpha, basis_norm2, basis_raw_coeffs);

        // Solve least squares: beta = (X^T X)^(-1) X^T y
        arma::mat XtX = X.t() * X;
        // Pseudo-inverse when singular; LAPACK-free in WASM builds.
        arma::mat XtX_inv = qe::linalg::inv_sym(XtX);

        arma::vec beta_x = XtX_inv * (X.t() * x);
        arma::vec beta_y = XtX_inv * (X.t() * y);

        arma::vec pred_x = X * beta_x;
        arma::vec pred_y = X * beta_y;

        arma::vec rx = x - pred_x;
        arma::vec ry = y - pred_y;

        // Hat matrix diagonal: h_ii = x_i^T (X^T X)^(-1) x_i
        arma::vec h(n);
        for (int i = 0; i < n; ++i) {
            arma::rowvec xi = X.row(i);
            h[i] = arma::as_scalar(xi * XtX_inv * xi.t());
            h[i] = std::min(0.9999, std::max(0.0, h[i]));
        }

        // 2D LOOCV error
        double cv_err = 0.0;
        for (int i = 0; i < n; ++i) {
            double denom = 1.0 - h[i];
            double rxi = rx[i] / denom;
            double ryi = ry[i] / denom;
            cv_err += (rxi * rxi + ryi * ryi);
        }

        // Bivariate AIC
        double rss = arma::dot(rx, rx) + arma::dot(ry, ry);
        double aic = static_cast<double>(n) * std::log(std::max(1e-12, rss / (2.0 * n))) + 2.0 * (2.0 * n_params);

        double score = (criterion == "aic") ? aic : cv_err;

        if (score < best_score || d == deg_start) {
            best_score = score;
            best_deg = d;
            best_fit.degree = d;
            best_fit.coeffs_x = basis_raw_coeffs * beta_x;
            best_fit.coeffs_y = basis_raw_coeffs * beta_y;
            best_fit.basis = basis;
            best_fit.basis_alpha = basis_alpha;
            best_fit.basis_norm2 = basis_norm2;
            best_fit.cv_error = cv_err;
            best_fit.aic = aic;
            best_fit.t = t_vec;
            best_fit.fitted = arma::join_rows(pred_x, pred_y);
        }
    }

    return best_fit;
}

/** @brief Compute trajectory velocities, speeds, accelerations, heading turn rates, and curvatures. */
inline TrajectoryDerivatives eval_trajectory_derivatives(
    const arma::vec& coeffs_x,
    const arma::vec& coeffs_y,
    const arma::vec& t_eval
) {
    size_t n = t_eval.n_elem;
    TrajectoryDerivatives d;
    d.t = t_eval;
    d.vx.set_size(n);
    d.vy.set_size(n);
    d.speed.set_size(n);
    d.ax.set_size(n);
    d.ay.set_size(n);
    d.heading_rate.set_size(n);
    d.curvature.set_size(n);

    for (size_t i = 0; i < n; ++i) {
        double t = t_eval[i];
        double vx = eval_poly_deriv1_scalar(coeffs_x, t);
        double vy = eval_poly_deriv1_scalar(coeffs_y, t);
        double ax = eval_poly_deriv2_scalar(coeffs_x, t);
        double ay = eval_poly_deriv2_scalar(coeffs_y, t);

        double spd = std::sqrt(vx * vx + vy * vy);
        double omega = vx * ay - vy * ax; // angular velocity / heading rate
        double spd3 = spd * spd * spd;
        double kappa = (spd3 > 1e-12) ? (omega / spd3) : 0.0;

        d.vx[i] = vx;
        d.vy[i] = vy;
        d.speed[i] = spd;
        d.ax[i] = ax;
        d.ay[i] = ay;
        d.heading_rate[i] = omega;
        d.curvature[i] = kappa;
    }

    return d;
}

} // namespace qe

#endif // LIBQE_TRAJECTORY_HPP

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <libqe/libqe.hpp>

using namespace Rcpp;

// =============================================================================
// Adjacency utilities
// =============================================================================

//' Number of upper-triangle pairs for n codes
//' @param n Number of codes
//' @return Integer: n*(n-1)/2
//' @examples
//' choose_two(4L)  # 6 code pairs
//' @export
// [[Rcpp::export]]
int choose_two(int n) {
    return qe::choose_two(n);
}

//' Number of undirected connection columns
//' @param n Number of codes
//' @param diagonal Include self-connections (A & A): n(n+1)/2 instead of choose(n, 2)
//' @examples
//' n_connections(4L)                   # 6: choose(4, 2)
//' n_connections(4L, diagonal = TRUE)  # 10: 4 * 5 / 2
//' @export
// [[Rcpp::export]]
int n_connections(int n, bool diagonal = false) {
    return qe::n_connections(n, diagonal);
}

//' Upper-triangle index pairs
//' @param len Number of codes (side length of square matrix)
//' @param row -1 = both rows, 0 = row indices only, 1 = col indices only
//' @param diagonal Include the diagonal (self-connections), in column-major
//'   order: (0,0), (0,1), (1,1), (0,2), ...
//' @examples
//' connection_indices(3L)                      # 2 x 3, 0-based (i, j) pairs
//' connection_indices(3L, row = 0L)            # row indices only
//' connection_indices(3L, diagonal = TRUE)     # includes (0,0), (1,1), (2,2)
//' @export
// [[Rcpp::export]]
arma::umat connection_indices(int len, int row = -1, bool diagonal = false) {
    return qe::connection_indices(len, row, diagonal);
}

//' Pairwise products as an upper-triangle vector
//' @param v Numeric matrix with one row of code values
//' @examples
//' code_connections(matrix(c(1, 0, 1), nrow = 1))  # A&B = 0, A&C = 1, B&C = 0
//' @export
// [[Rcpp::export]]
arma::rowvec code_connections(arma::mat v) {
    return qe::code_connections(v);
}

//' Fold a directed (n*n) vector into an undirected upper-triangle vector
//' @param v Numeric vector of length n*n
//' @param diagonal Keep the diagonal (self-connections, each once): A & A, A & B, B & B, ...
//' @examples
//' m <- matrix(c(0, 1, 2,
//'               3, 0, 4,
//'               5, 6, 0), nrow = 3, byrow = TRUE)
//' fold_directed_network(as.vector(m))                   # A&B, A&C, B&C
//' fold_directed_network(as.vector(m), diagonal = TRUE)
//' @export
// [[Rcpp::export]]
arma::rowvec fold_directed_network(arma::vec v, bool diagonal = false) {
    return qe::fold_directed_network(v, diagonal);
}

//' Flatten an adjacency matrix to a connection vector
//' @param x Numeric matrix
//' @param full TRUE = full n*n (directed); FALSE = upper triangle (undirected)
//' @param diagonal With full = FALSE: include the diagonal (A & A, A & B, B & B, ...)
//' @examples
//' m <- matrix(c(0, 1, 2,
//'               1, 0, 3,
//'               2, 3, 0), nrow = 3)
//' network_to_vector(m, full = FALSE)  # upper triangle: 1 2 3
//' network_to_vector(m, full = TRUE)   # all 9 entries
//' @export
// [[Rcpp::export]]
arma::rowvec network_to_vector(arma::mat x, bool full = true, bool diagonal = false) {
    return qe::network_to_vector(x, full, diagonal);
}

//' Code-name pairs for upper-triangle positions ("A & B")
//' @param v Character vector of code names
//' @param diagonal Include self-connections: "A & A", "A & B", "B & B", ...
//' @examples
//' connection_names(c("A", "B", "C"))
//' connection_names(c("A", "B", "C"), diagonal = TRUE)
//' @export
// [[Rcpp::export]]
std::vector<std::string> connection_names(std::vector<std::string> v, bool diagonal = false) {
    return qe::connection_names(v, diagonal);
}

// =============================================================================
// Normalization
// =============================================================================

//' Row-wise L2 (sphere) normalization
//' @param m Numeric matrix
//' @examples
//' m <- matrix(c(3, 4, 0,
//'               0, 0, 0,
//'               1, 1, 1), nrow = 3, byrow = TRUE)
//' normalize_networks(m)  # rows have unit length; the zero row is unchanged
//' @export
// [[Rcpp::export]]
arma::mat normalize_networks(arma::mat m) {
    return qe::normalize_networks(m);
}

//' Max-norm scaling (divide all rows by the largest row L2 norm)
//' @param m Numeric matrix
//' @examples
//' m <- matrix(c(3, 4,
//'               1, 0), nrow = 2, byrow = TRUE)
//' scale_networks(m)  # divided by 5, the largest row norm
//' @export
// [[Rcpp::export]]
arma::mat scale_networks(arma::mat m) {
    return qe::scale_networks(m);
}

// =============================================================================
// Modeling
// =============================================================================

//' Center points (subtract column means)
//' @param values Numeric matrix
//' @examples
//' pts <- matrix(c(1, 2, 3, 10, 20, 30), ncol = 2)
//' center_points(pts)  # column means are now 0
//' @export
// [[Rcpp::export]]
arma::mat center_points(arma::mat values) {
    return qe::center_points(values);
}

//' Confidence interval for the mean of a group of ENA unit points
//'
//' For each dimension computes a t-based CI:
//' \code{mean +/- t(alpha/2, n-1) * (SD / sqrt(n))} where \code{alpha = 1 - conf_level}.
//'
//' @param points     Numeric matrix (units x dims), one row per unit in the group
//' @param conf_level Confidence level (default 0.95)
//' @return Numeric matrix (n_dims x 3): columns are [mean, ci_lower, ci_upper].
//'   When \code{nrow(points) == 1} the CI bounds are \code{-Inf} and \code{Inf}.
//' @examples
//' pts <- matrix(c(1, 2, 3, 4, 10, 12, 11, 13), ncol = 2)
//' mean_ci(pts)  # rows = dimensions; columns = mean, lower, upper
//' t.test(pts[, 1])$conf.int  # same bounds as row 1
//' @export
// [[Rcpp::export]]
arma::mat mean_ci(arma::mat points, double conf_level = 0.95) {
    return qe::mean_ci(points, conf_level);
}

//' Outlier interval based on IQR (Tukey fence) for a group of ENA unit points
//'
//' For each dimension: \code{lower = -iqr_factor * IQR}, \code{upper = +iqr_factor * IQR}.
//' Equivalent to rENA's:
//' \preformatted{
//'   oi <- c(IQR(pts[,1]), ...) * iqr_factor
//'   matrix(rep(oi, 2), ncol = n_dims, byrow = TRUE) * c(-1, 1)
//' }
//' IQR uses R's default \code{type = 7} quantile.
//'
//' @param points     Numeric matrix (units x dims), one row per unit in the group
//' @param iqr_factor Multiplier applied to the IQR (default 1.5, the Tukey fence)
//' @return Numeric matrix (n_dims x 2): columns are [lower, upper].
//' @examples
//' pts <- matrix(c(1, 2, 3, 4, 10, 12, 11, 13), ncol = 2)
//' outlier_ci(pts)  # rows = dimensions; columns = -1.5 * IQR, +1.5 * IQR
//' @export
// [[Rcpp::export]]
arma::mat outlier_ci(arma::mat points, double iqr_factor = 1.5) {
    return qe::outlier_ci(points, iqr_factor);
}

//' Two-group comparison statistics
//'
//' Computes per-dimension parametric and non-parametric statistics for two
//' groups of ENA unit points.  Equivalent to \code{group.stats()} in rENA-api,
//' moved to the C++ layer so all language bindings share a single implementation.
//'
//' \strong{Parametric (Welch two-sample t-test):}
//' \itemize{
//'   \item \code{t}        - Welch t-statistics
//'   \item \code{parameter} - Welch-Satterthwaite degrees of freedom
//'   \item \code{pvalue}   - two-tailed p-values
//'   \item \code{effect}   - Cohen's d (pooled SD: \code{(mean1 - mean2) / pooled_sd})
//'   \item \code{mean}     - 2 x n_dims matrix, row 1 = group1, row 2 = group2
//'   \item \code{std.dev}  - 2 x n_dims matrix of sample standard deviations
//' }
//'
//' \strong{Non-parametric (Wilcoxon rank-sum):}
//' \itemize{
//'   \item \code{U}        - rank-sum U for group 1 (= R's \code{wilcox.test} W)
//'   \item \code{pvalue}   - two-tailed p-values (normal approx, tie + continuity correction)
//'   \item \code{effect}   - rank-biserial: \code{1 - 2 * U / (n1 * n2)}
//'   \item \code{median}   - 2 x n_dims matrix of medians
//' }
//'
//' Entries are \code{NA} when a statistic is undefined (e.g. n < 2 for
//' parametric tests).
//'
//' @param g1 Numeric matrix (n1 x n_dims), unit points for group 1
//' @param g2 Numeric matrix (n2 x n_dims), unit points for group 2
//' @return Named list with \code{N}, \code{parametric}, and \code{nonparametric}
//'   sub-lists, matching the structure of rENA-api's \code{group.stats()}.
//' @examples
//' g1 <- matrix(c(1, 2, 3, 4, 10, 11, 12, 13), ncol = 2)
//' g2 <- matrix(c(2, 3, 5, 6, 9, 9, 10, 12), ncol = 2)
//' s <- group_stats(g1, g2)
//' s$N
//' s$parametric$t
//' s$nonparametric$U
//' @export
// [[Rcpp::export]]
List group_stats(arma::mat g1, arma::mat g2) {
    qe::GroupStats s = qe::group_stats(g1, g2);
    return List::create(
        _("N") = IntegerVector::create(s.n1, s.n2),
        _("parametric") = List::create(
            _("t")         = s.t,
            _("parameter") = s.df,
            _("pvalue")    = s.pvalue_t,
            _("effect")    = s.cohens_d,
            _("mean")      = s.means,
            _("std.dev")   = s.sds
        ),
        _("nonparametric") = List::create(
            _("U")      = s.U,
            _("pvalue") = s.pvalue_u,
            _("effect") = s.effect_r,
            _("median") = s.medians
        )
    );
}

// =============================================================================
// Door Temporal Pooling & Smoothing
// =============================================================================

//' Door lookback window pooling for a single block
//' @param block Numeric matrix of connection counts
//' @param lookback_size Window lookback size
//' @param aggregate_mean Logical: TRUE = mean, FALSE = sum
//' @param weighting_linear Logical: TRUE = linear weights, FALSE = equal
//' @param segment_ids Integer vector of segment IDs (0-indexed or 1-indexed)
//' @examples
//' block <- matrix(c(1, 0,
//'                   2, 1,
//'                   3, 0,
//'                   4, 2), nrow = 4, byrow = TRUE)
//' door_lookback_block(block, lookback_size = 2L)                        # sums
//' door_lookback_block(block, lookback_size = 2L, aggregate_mean = TRUE) # means
//' @export
// [[Rcpp::export]]
arma::mat door_lookback_block(
    arma::mat block,
    int lookback_size = 20,
    bool aggregate_mean = false,
    bool weighting_linear = false,
    IntegerVector segment_ids = IntegerVector()
) {
    std::vector<int> segs(segment_ids.begin(), segment_ids.end());
    return qe::lookback_block(block, lookback_size, aggregate_mean, weighting_linear, segs);
}

//' Door EMA smoothing for a single block
//' @param block Numeric matrix of connection counts
//' @param alpha Smoothing factor in (0, 1]
//' @param segment_ids Integer vector of segment IDs
//' @examples
//' block <- matrix(c(10, 0,
//'                    0, 10,
//'                    0,  0), nrow = 3, byrow = TRUE)
//' door_ema_block(block, alpha = 0.5)
//' @export
// [[Rcpp::export]]
arma::mat door_ema_block(
    arma::mat block,
    double alpha = 0.1,
    IntegerVector segment_ids = IntegerVector()
) {
    std::vector<int> segs(segment_ids.begin(), segment_ids.end());
    return qe::ema_block(block, alpha, segs);
}

//' Door lookback window pooling across units
//' @param conn_counts Numeric matrix of connection counts
//' @param unit_row_indices List of integer vectors of row indices (0-indexed or 1-indexed R converted)
//' @param lookback_sizes Integer vector of lookback sizes per unit
//' @param aggregate_mean Logical: TRUE = mean, FALSE = sum
//' @param weighting_linear Logical: TRUE = linear weights, FALSE = equal
//' @param segment_ids Optional integer vector of segment IDs
//' @examples
//' counts <- matrix(c(1, 0,
//'                    2, 1,
//'                    5, 5,
//'                    6, 5), nrow = 4, byrow = TRUE)
//' # two units, rows given 0-based; each unit pools within its own rows
//' door_lookback(counts, list(c(0L, 1L), c(2L, 3L)), lookback_sizes = c(2L, 2L))
//' @export
// [[Rcpp::export]]
arma::mat door_lookback(
    arma::mat conn_counts,
    List unit_row_indices,
    IntegerVector lookback_sizes,
    bool aggregate_mean = false,
    bool weighting_linear = false,
    IntegerVector segment_ids = IntegerVector()
) {
    std::vector<std::vector<int>> unit_rows;
    unit_rows.reserve(unit_row_indices.size());
    for (int i = 0; i < unit_row_indices.size(); ++i) {
        IntegerVector iv = unit_row_indices[i];
        unit_rows.push_back(std::vector<int>(iv.begin(), iv.end()));
    }
    std::vector<int> lbs(lookback_sizes.begin(), lookback_sizes.end());
    std::vector<int> segs(segment_ids.begin(), segment_ids.end());

    return qe::apply_door_lookback(conn_counts, unit_rows, lbs, aggregate_mean, weighting_linear, segs);
}

//' Door EMA smoothing across units
//' @param conn_counts Numeric matrix of connection counts
//' @param unit_row_indices List of integer vectors of row indices
//' @param alpha Smoothing factor
//' @param segment_ids Optional integer vector of segment IDs
//' @examples
//' counts <- matrix(c(10, 0,
//'                     0, 10,
//'                    10, 0,
//'                     0, 10), nrow = 4, byrow = TRUE)
//' door_ema(counts, list(c(0L, 1L), c(2L, 3L)), alpha = 0.5)
//' @export
// [[Rcpp::export]]
arma::mat door_ema(
    arma::mat conn_counts,
    List unit_row_indices,
    double alpha = 0.1,
    IntegerVector segment_ids = IntegerVector()
) {
    std::vector<std::vector<int>> unit_rows;
    unit_rows.reserve(unit_row_indices.size());
    for (int i = 0; i < unit_row_indices.size(); ++i) {
        IntegerVector iv = unit_row_indices[i];
        unit_rows.push_back(std::vector<int>(iv.begin(), iv.end()));
    }
    std::vector<int> segs(segment_ids.begin(), segment_ids.end());

    return qe::apply_door_ema(conn_counts, unit_rows, alpha, segs);
}

// =============================================================================
// Trajectory Curve Fitting & Differential Geometry
// =============================================================================

//' Fit 2D parametric polynomial curve with LOOCV-2D or AIC degree selection
//' @param points Numeric matrix (n x 2) of trajectory coordinates
//' @param t Numeric vector of time values (length n)
//' @param max_degree Maximum polynomial degree (default 3)
//' @param fixed_degree Exact degree to use when >= 1
//' @param criterion Criterion for degree selection: "loocv" or "aic"
//' @param basis Polynomial basis used for fitting: "orthogonal" or "raw"
//' @examples
//' t <- seq(0, 1, length.out = 6)
//' pts <- cbind(1 + 2 * t, 4 - 3 * t)  # a straight line
//' fit <- fit_trajectory_poly(pts, t = t, fixed_degree = 1L, basis = "raw")
//' fit$degree
//' fit$coeffs_x  # intercept 1, slope 2
//' fit$coeffs_y  # intercept 4, slope -3
//' @export
// [[Rcpp::export]]
List fit_trajectory_poly(
    arma::mat points,
    NumericVector t = NumericVector(),
    int max_degree = 3,
    int fixed_degree = 0,
    std::string criterion = "loocv",
    std::string basis = "orthogonal"
) {
    arma::vec t_vec;
    if (t.size() > 0) {
        t_vec = arma::vec(t.begin(), t.size());
    }
    qe::PolyCurveFit fit = qe::fit_poly_loocv2d(points, t_vec, max_degree, fixed_degree, criterion, basis);

    return List::create(
        Named("degree")        = fit.degree,
        Named("coeffs_x")      = fit.coeffs_x,
        Named("coeffs_y")      = fit.coeffs_y,
        Named("basis")         = fit.basis,
        Named("basis_alpha")   = fit.basis_alpha,
        Named("basis_norm2")   = fit.basis_norm2,
        Named("cv_error")      = fit.cv_error,
        Named("aic")           = fit.aic,
        Named("t")             = fit.t,
        Named("fitted_points") = fit.fitted
    );
}

//' Evaluate 2D parametric polynomial curve at specified time points
//' @param coeffs_x Polynomial coefficients for X
//' @param coeffs_y Polynomial coefficients for Y
//' @param t_eval Numeric vector of time evaluation points
//' @examples
//' # x(t) = 1 + 2 t, y(t) = 4 - 3 t (raw polynomial coefficients)
//' eval_trajectory_curve(c(1, 2), c(4, -3), c(0, 0.5, 1))
//' @export
// [[Rcpp::export]]
arma::mat eval_trajectory_curve(
    arma::vec coeffs_x,
    arma::vec coeffs_y,
    arma::vec t_eval
) {
    return qe::eval_poly_curve(coeffs_x, coeffs_y, t_eval);
}

//' Evaluate trajectory differential geometry (velocities, speed, curvature, turns)
//' @param coeffs_x Polynomial coefficients for X
//' @param coeffs_y Polynomial coefficients for Y
//' @param t_eval Numeric vector of time evaluation points
//' @examples
//' # x(t) = t, y(t) = t^2
//' d <- eval_trajectory_derivatives(c(0, 1), c(0, 0, 1), c(0, 0.5, 1))
//' d$speed
//' d$curvature
//' @export
// [[Rcpp::export]]
List eval_trajectory_derivatives(
    arma::vec coeffs_x,
    arma::vec coeffs_y,
    arma::vec t_eval
) {
    qe::TrajectoryDerivatives d = qe::eval_trajectory_derivatives(coeffs_x, coeffs_y, t_eval);
    return List::create(
        Named("t")            = d.t,
        Named("vx")           = d.vx,
        Named("vy")           = d.vy,
        Named("speed")        = d.speed,
        Named("ax")           = d.ax,
        Named("ay")           = d.ay,
        Named("heading_rate") = d.heading_rate,
        Named("curvature")    = d.curvature
    );
}

// =============================================================================
// Trajectory Distance & Following Dynamics
// =============================================================================

//' Integrated Euclidean distance between two trajectory polynomial curves
//' @param coeffs_ax Polynomial coeffs for Curve A X
//' @param coeffs_ay Polynomial coeffs for Curve A Y
//' @param coeffs_bx Polynomial coeffs for Curve B X
//' @param coeffs_by Polynomial coeffs for Curve B Y
//' @param t_start Integration start (default 0.0)
//' @param t_end Integration end (default 1.0)
//' @examples
//' cx <- c(0, 1); cy <- c(0, 2)
//' integrated_trajectory_distance(cx, cy, cx, cy)              # 0: same curve
//' integrated_trajectory_distance(cx, cy, cx + c(3, 0), cy + c(4, 0))  # 5: parallel offset
//' @export
// [[Rcpp::export]]
double integrated_trajectory_distance(
    arma::vec coeffs_ax,
    arma::vec coeffs_ay,
    arma::vec coeffs_bx,
    arma::vec coeffs_by,
    double t_start = 0.0,
    double t_end = 1.0
) {
    return qe::integrated_curve_distance(coeffs_ax, coeffs_ay, coeffs_bx, coeffs_by, t_start, t_end);
}

//' Lagged curve distance between follower and leader
//' @param coeffs_fol_x Polynomial coeffs for Follower X
//' @param coeffs_fol_y Polynomial coeffs for Follower Y
//' @param coeffs_ldr_x Polynomial coeffs for Leader X
//' @param coeffs_ldr_y Polynomial coeffs for Leader Y
//' @param lag Time lag in [0, 1)
//' @examples
//' cx <- c(0, 1); cy <- c(0, 2)
//' lagged_trajectory_distance(cx, cy, cx, cy, lag = 0)
//' lagged_trajectory_distance(cx, cy, cx, cy, lag = 0.25)
//' @export
// [[Rcpp::export]]
double lagged_trajectory_distance(
    arma::vec coeffs_fol_x,
    arma::vec coeffs_fol_y,
    arma::vec coeffs_ldr_x,
    arma::vec coeffs_ldr_y,
    double lag = 0.0
) {
    return qe::lagged_curve_distance(coeffs_fol_x, coeffs_fol_y, coeffs_ldr_x, coeffs_ldr_y, lag);
}

//' Pairwise integrated trajectory distance matrix
//' @param all_coeffs_x List of X coefficient vectors
//' @param all_coeffs_y List of Y coefficient vectors
//' @examples
//' cx <- list(c(0, 1), c(3, 1), c(0, 1))
//' cy <- list(c(0, 2), c(4, 2), c(0, 2))
//' pairwise_trajectory_distance(cx, cy)  # 3 x 3; curves 1 and 3 coincide
//' @export
// [[Rcpp::export]]
arma::mat pairwise_trajectory_distance(
    List all_coeffs_x,
    List all_coeffs_y
) {
    std::vector<arma::vec> cx;
    std::vector<arma::vec> cy;
    cx.reserve(all_coeffs_x.size());
    cy.reserve(all_coeffs_y.size());
    for (int i = 0; i < all_coeffs_x.size(); ++i) {
        NumericVector nx = all_coeffs_x[i];
        NumericVector ny = all_coeffs_y[i];
        cx.push_back(arma::vec(nx.begin(), nx.size()));
        cy.push_back(arma::vec(ny.begin(), ny.size()));
    }
    return qe::pairwise_trajectory_distance_matrix(cx, cy);
}

//' Signed turn-lag distance between two agents at specific lag delta
//' @param pts_a Points matrix for agent A (n_a x 2)
//' @param pts_b Points matrix for agent B (n_b x 2)
//' @param times_a Turn times for agent A
//' @param times_b Turn times for agent B
//' @param delta Signed lag in turns
//' @examples
//' times_b <- c(1, 2, 3, 4, 5)
//' pts_b <- cbind(1:5, 1:5)
//' # agent A repeats B's positions two turns later
//' times_a <- times_b + 2
//' pts_a <- pts_b
//' signed_turn_lag(pts_a, pts_b, times_a, times_b, delta = 2L)  # mean distance 0
//' signed_turn_lag(pts_a, pts_b, times_a, times_b, delta = 0L)
//' @export
// [[Rcpp::export]]
List signed_turn_lag(
    arma::mat pts_a,
    arma::mat pts_b,
    arma::vec times_a,
    arma::vec times_b,
    int delta
) {
    auto [mean_d, count] = qe::signed_turn_lag_distance(pts_a, pts_b, times_a, times_b, delta);
    return List::create(
        Named("mean_distance") = mean_d,
        Named("valid_count")   = count
    );
}

//' Sweep signed turn lags to find optimal leader-follower lag
//' @param pts_a Points matrix for agent A (n_a x 2)
//' @param pts_b Points matrix for agent B (n_b x 2)
//' @param times_a Turn times for agent A
//' @param times_b Turn times for agent B
//' @param max_lag Maximum turn lag to evaluate (default 15)
//' @examples
//' times_b <- c(1, 2, 3, 4, 5)
//' pts_b <- cbind(1:5, 1:5)
//' times_a <- times_b + 2  # A follows B by two turns
//' pts_a <- pts_b
//' res <- sweep_signed_turn_lags(pts_a, pts_b, times_a, times_b, max_lag = 4L)
//' res$best_lag  # 2
//' res$mean_distances
//' @export
// [[Rcpp::export]]
List sweep_signed_turn_lags(
    arma::mat pts_a,
    arma::mat pts_b,
    arma::vec times_a,
    arma::vec times_b,
    int max_lag = 15
) {
    qe::SignedTurnLagResult res = qe::best_signed_turn_lag(pts_a, pts_b, times_a, times_b, max_lag);
    return List::create(
        Named("best_lag")           = res.best_lag,
        Named("min_mean_distance")  = res.min_mean_distance,
        Named("lags")               = res.lags,
        Named("mean_distances")     = res.mean_distances,
        Named("valid_counts")       = res.valid_counts
    );
}

//' Distance-distance matrix correlation for stability analysis
//' @param X Numeric matrix
//' @param Y Numeric matrix
//' @examples
//' X <- matrix(c(0, 1, 0, 1, 0, 0, 1, 1), ncol = 2)
//' Y <- 2 * X  # a uniform scaling keeps the distance structure
//' dist_dist_correlation(X, Y)  # 1
//' @export
// [[Rcpp::export]]
double dist_dist_correlation(
    arma::mat X,
    arma::mat Y
) {
    return qe::dist_dist_correlation(X, Y);
}

1		#' @useDynLib libqe, .registration = TRUE
2		#' @importFrom Rcpp evalCpp
3		NULL
4
5		.onLoad <- function(libname, pkgname) {
6		# Nothing to initialise: libqe is a header + Rcpp library.
7	!	invisible(NULL)
8		}

1		/**
2		* @file adjacency.hpp
3		* @brief Combinatorics and vector/matrix utilities for ENA connection networks.
4		*
5		* Pure C++ / Armadillo — no Rcpp dependency.
6		* Include @c <RcppArmadillo.h> before this header when building inside an R package.
7		*/
8		#ifndef LIBQE_ADJACENCY_HPP
9		#define LIBQE_ADJACENCY_HPP
10
11		#include <armadillo>
12		#include <string>
13		#include <vector>
14		#include <cmath>
15
16		namespace qe {
17
18		// ---------------------------------------------------------------------------
19		/// @name Combinatorics
20		/// @{
21		// ---------------------------------------------------------------------------
22
23		/**
24		* @brief Number of unordered pairs from @p n items (n choose 2).
25		*
26		* @param n Number of codes/nodes.
27		* @returns @c (n * (n - 1)) / 2.
28		*/
29	15x	inline int choose_two(int n) {
30	15x	return (n * (n - 1)) / 2;
31		}
32
33		/**
34		* @brief Number of undirected connection columns for @p n codes.
35		*
36		* @param n Number of codes/nodes.
37		* @param diagonal Include self-connections (the diagonal, @c A & A).
38		* @returns @c choose_two(n), or @c n(n+1)/2 with @p diagonal.
39		*/
40	16x	inline int n_connections(int n, bool diagonal = false) {
41	16x	return diagonal ? (n * (n + 1)) / 2 : choose_two(n);
42		}
43
44		/**
45		* @brief Upper-triangle index pairs for a symmetric matrix of side @p len.
46		*
47		* Column-major order: (0,1), (0,2), (1,2), ... -- or, with @p diagonal,
48		* (0,0), (0,1), (1,1), (0,2), (1,2), (2,2), ... (Armadillo's
49		* @c trimatu_ind order), so a self-connection's row and column index are
50		* equal.
51		*
52		* @param len Side length of the square matrix (number of codes).
53		* @param row Controls what is returned:
54		* - @c -1 (default): 2 × k matrix of [row_idx; col_idx].
55		* - @c 0: row indices only (1 × k).
56		* - @c 1: column indices only (1 × k).
57		* @param diagonal Include the diagonal (self-connections). Default @c false.
58		* @returns An @c arma::umat of shape 2×k or 1×k depending on @p row, with
59		* @c k = n_connections(len, diagonal).
60		* @note Equivalent to @c triIndices() in rENA/ena.cpp and tma/code.cpp
61		* (without @p diagonal).
62		*/
63	10x	inline arma::umat connection_indices(int len, int row = -1, bool diagonal = false) {
64	10x	int vS = n_connections(len, diagonal);
65	10x	int s = 0;
66
67	10x	arma::umat vR(2, vS, arma::fill::zeros);
68	10x	arma::umat vRone(1, vS, arma::fill::zeros);
69
70	58x	for (int i = 1; i <= len; i++) {
71	163x	for (int j = 0; j < (diagonal ? i : i - 1); j++) {
72	115x	vR(0, s) = j;
73	115x	vR(1, s) = i - 1;
74	121x	if (row == 0) vRone[s] = j;
75	115x	else if (row == 1) vRone[s] = i - 1;
76	115x	s++;
77		}
78		}
79	20x	return (row == -1) ? vR : vRone;
80		}
81
82		/// @}
83
84		// ---------------------------------------------------------------------------
85		/// @name Vector ↔ upper-triangle conversions
86		/// @{
87		// ---------------------------------------------------------------------------
88
89		/**
90		* @brief Compute pairwise products and return as a flat upper-triangle vector.
91		*
92		* For each pair @c (j, i) with @c j < i, outputs @c v[j] * v[i].
93		* The result has length @c choose_two(v.size()).
94		*
95		* @param v A 1 × p row vector (or p-element matrix).
96		* @returns A row vector of length @c choose_two(p).
97		* @note Equivalent to @c vector_to_ut() in rENA/ena.cpp.
98		*/
99	2x	inline arma::rowvec code_connections(arma::mat v) {
100	2x	int vL = v.size();
101	2x	int vS = choose_two(vL);
102	2x	arma::rowvec out(vS, arma::fill::zeros);
103	2x	int s = 0;
104	7x	for (int i = 2; i <= vL; i++) {
105	14x	for (int j = 0; j < i - 1; j++) {
106	27x	out[s] = v[j] * v[i - 1];
107	9x	s++;
108		}
109		}
110	2x	return out;
111		}
112
113		/**
114		* @brief Fold a directed n² vector into an undirected upper-triangle vector.
115		*
116		* Symmetric elements are summed: A→B + B→A for every pair. With
117		* @p diagonal the self-connections (A→A) are kept, each once, in the
118		* @c connection_indices(n, -1, true) order (A & A, A & B, B & B, ...).
119		*
120		* @param v Flat directed vector of length n², column-major
121		* (n = sqrt(v.size())).
122		* @param diagonal Keep the diagonal. Default @c false.
123		* @returns Upper-triangle row vector of length @c n_connections(n, diagonal).
124		* @note Equivalent to @c vector_to_summed_uppertri() in tma/code.cpp
125		* (without @p diagonal).
126		*/
127	7x	inline arma::rowvec fold_directed_network(arma::vec v, bool diagonal = false) {
128	7x	int n = static_cast<int>(std::round(std::sqrt(static_cast<double>(v.size()))));
129
130	7x	arma::mat m(v);
131	7x	m.reshape(n, n);
132	14x	arma::mat folded = m + m.t();
133	9x	if (diagonal) folded.diag() = m.diag(); // a self-connection is not a pair
134
135	7x	arma::colvec flat = folded.as_col();
136	7x	arma::uvec up_inds = arma::trimatu_ind(arma::size(folded), diagonal ? 0 : 1);
137	14x	return arma::conv_to<arma::rowvec>::from(flat.elem(up_inds));
138		}
139
140		/**
141		* @brief Flatten an adjacency matrix to a vector.
142		*
143		* @param x Square adjacency matrix (n × n).
144		* @param full If @c true (default), returns the full n² vector (directed).
145		* If @c false, returns the upper-triangle only (undirected).
146		* @param diagonal With @p full @c false: include the diagonal, in the
147		* @c connection_indices(n, -1, true) order. Default @c false.
148		* @returns Row vector of length n² or @c n_connections(n, diagonal).
149		* @note Equivalent to @c adjacency_matrix_to_vector() in tma/code.cpp
150		* (without @p diagonal).
151		*/
152	3x	inline arma::rowvec network_to_vector(arma::mat x, bool full = true, bool diagonal = false) {
153	4x	if (full) return arma::vectorise(x).t();
154	2x	if (diagonal) {
155	1x	arma::uvec up_inds = arma::trimatu_ind(arma::size(x), 0);
156	2x	return arma::vectorise(x(up_inds)).t();
157		}
158	3x	arma::mat combined = arma::trimatu(x, 1) + arma::trimatl(x, -1);
159	1x	arma::uvec up_inds = arma::trimatu_ind(arma::size(combined), 1);
160	2x	return arma::vectorise(combined(up_inds)).t();
161		}
162
163		/// @}
164
165		// ---------------------------------------------------------------------------
166		/// @name String utilities
167		/// @{
168		// ---------------------------------------------------------------------------
169
170		/**
171		* @brief Generate "A & B" pair names for every upper-triangle position.
172		*
173		* For code names @c {"A", "B", "C"} returns @c {"A & B", "A & C", "B & C"};
174		* with @p diagonal, @c {"A & A", "A & B", "B & B", "A & C", "B & C", "C & C"}
175		* (the @c connection_indices(p, -1, true) order).
176		*
177		* @param v Ordered list of code names (length p).
178		* @param diagonal Include self-connections. Default @c false.
179		* @returns Vector of @c n_connections(p, diagonal) label strings.
180		* @note Equivalent to @c svector_to_ut() in rENA/ena.cpp (without @p diagonal).
181		*/
182	4x	inline std::vector<std::string> connection_names(std::vector<std::string> v,
183		bool diagonal = false) {
184	4x	int vL = v.size();
185	4x	int vS = n_connections(vL, diagonal);
186	4x	std::vector<std::string> out(vS);
187	4x	int s = 0;
188	17x	for (int i = 1; i <= vL; i++) {
189	31x	for (int j = 0; j < (diagonal ? i : i - 1); j++) {
190	18x	out[s] = v[j] + " & " + v[i - 1];
191	18x	s++;
192		}
193		}
194	4x	return out;
195		}
196
197		/// @}
198
199		} // namespace qe
200
201		#endif // LIBQE_ADJACENCY_HPP

1		/** @file door.hpp
2		* @brief Door temporal pooling and smoothing functions for trajectory modeling.
3		*
4		* Provides sliding lookback window pooling (with uniform or linear weighting,
5		* sum or mean aggregation, and optional segment boundary resets) as well as
6		* Exponential Moving Average (EMA) smoothing over consecutive time steps.
7		*/
8		#ifndef LIBQE_DOOR_HPP
9		#define LIBQE_DOOR_HPP
10
11		#include <armadillo>
12		#include <vector>
13		#include <string>
14		#include <cmath>
15		#include <algorithm>
16
17		namespace qe {
18
19		/** @brief Apply a sliding lookback window pool over a single unit's connection matrix.
20		*
21		* For each row @c i, pools rows in the window [@c max(0, i - lookback_size + 1), @c i].
22		* If @p segment_ids is non-empty, the window does not cross rows with a different segment ID.
23		*
24		* @param[in] block Matrix of connection counts for one unit (n_rows × n_features).
25		* @param[in] lookback_size Number of rows in the lookback window (>= 1).
26		* @param[in] aggregate_mean When @c true, averages rows; when @c false, sums rows.
27		* @param[in] weighting_linear When @c true, applies linearly increasing weights
28		* (oldest row in window has weight 1, current row has weight k).
29		* When @c false, all rows in window have equal weight 1.
30		* @param[in] segment_ids Optional vector of segment/activity IDs (length n_rows).
31		* Window pooling resets at segment transitions.
32		*
33		* @returns A matrix of shape n_rows × n_features containing door-smoothed values.
34		*/
35	5x	inline arma::mat lookback_block(
36		const arma::mat& block,
37		int lookback_size = 20,
38		bool aggregate_mean = false,
39		bool weighting_linear = false,
40		const std::vector<int>& segment_ids = {}
41		) {
42	5x	int n_rows = static_cast<int>(block.n_rows);
43	5x	int n_cols = static_cast<int>(block.n_cols);
44	5x	arma::mat out(n_rows, n_cols, arma::fill::zeros);
45
46		if (n_rows == 0 \|\| n_cols == 0) return out;
47	5x	int k_max = std::max(1, lookback_size);
48	5x	bool has_segments = (!segment_ids.empty() && static_cast<int>(segment_ids.size()) >= n_rows);
49
50	25x	for (int i = 0; i < n_rows; ++i) {
51	20x	int start_idx = std::max(0, i - k_max + 1);
52
53	20x	if (has_segments) {
54	4x	int cur_seg = segment_ids[i];
55	10x	for (int s = i; s >= start_idx; --s) {
56	8x	if (segment_ids[s] != cur_seg) {
57	2x	start_idx = s + 1;
58	2x	break;
59		}
60		}
61		}
62
63	20x	int win_len = i - start_idx + 1;
64		if (win_len <= 0) continue;
65
66	20x	arma::rowvec acc(n_cols, arma::fill::zeros);
67	20x	arma::rowvec denom(n_cols, arma::fill::zeros);
68
69	56x	for (int j = start_idx; j <= i; ++j) {
70	36x	double w = 1.0;
71	36x	if (weighting_linear) {
72		// Linear weight: 1 for earliest in window, win_len for current row i
73	!	w = static_cast<double>(j - start_idx + 1);
74		}
75	108x	for (int col = 0; col < n_cols; ++col) {
76	72x	double value = block(j, col);
77	72x	if (!std::isnan(value)) {
78	66x	acc[col] += w * value;
79	132x	denom[col] += w;
80		}
81		}
82		}
83
84	20x	if (aggregate_mean) {
85	24x	for (int col = 0; col < n_cols; ++col) {
86		acc[col] = (denom[col] > 0.0) ? (acc[col] / denom[col]) : arma::datum::nan;
87		}
88		}
89
90	40x	out.row(i) = acc;
91		}
92
93	5x	return out;
94		}
95
96		/** @brief Apply Exponential Moving Average (EMA) smoothing over a single unit's connection matrix.
97		*
98		* Computes S_0 = Y_0, S_t = alpha * Y_t + (1 - alpha) * S_{t-1}.
99		* If @p segment_ids is non-empty, S_t resets to Y_t when transitioning to a new segment.
100		*
101		* @param[in] block Matrix of connection counts for one unit (n_rows × n_features).
102		* @param[in] alpha Smoothing parameter alpha in (0, 1].
103		* @param[in] segment_ids Optional vector of segment/activity IDs (length n_rows).
104		*
105		* @returns A matrix of shape n_rows × n_features containing EMA-smoothed values.
106		*/
107	2x	inline arma::mat ema_block(
108		const arma::mat& block,
109		double alpha = 0.1,
110		const std::vector<int>& segment_ids = {}
111		) {
112	2x	int n_rows = static_cast<int>(block.n_rows);
113	2x	int n_cols = static_cast<int>(block.n_cols);
114	2x	arma::mat out(n_rows, n_cols, arma::fill::zeros);
115
116		if (n_rows == 0 \|\| n_cols == 0) return out;
117	2x	double a = std::max(1e-6, std::min(1.0, alpha));
118		bool has_segments = (!segment_ids.empty() && static_cast<int>(segment_ids.size()) >= n_rows);
119
120	2x	out.row(0) = block.row(0);
121
122	6x	for (int i = 1; i < n_rows; ++i) {
123		bool reset = has_segments && (segment_ids[i] != segment_ids[i - 1]);
124	4x	if (reset) {
125	!	out.row(i) = block.row(i);
126		} else {
127	12x	for (int col = 0; col < n_cols; ++col) {
128	8x	double cur = block(i, col);
129	8x	if (std::isnan(cur)) {
130	4x	cur = out(i - 1, col);
131		}
132	24x	out(i, col) = a * cur + (1.0 - a) * out(i - 1, col);
133		}
134		}
135		}
136
137	2x	return out;
138		}
139
140		/** @brief Apply lookback door pooling to an entire connection matrix across multiple units.
141		*
142		* @param[in] conn_counts Complete connection matrix (all units).
143		* @param[in] unit_row_indices Vector of row-index vectors, one per unit (0-indexed).
144		* @param[in] lookback_sizes Vector of lookback sizes, one per unit.
145		* @param[in] aggregate_mean When @c true, compute mean; when @c false, compute sum.
146		* @param[in] weighting_linear When @c true, use linear weighting.
147		* @param[in] segment_ids Optional vector of segment IDs for every row in conn_counts.
148		*
149		* @returns Smoothed connection matrix with identical dimensions to @p conn_counts.
150		*/
151	!	inline arma::mat apply_door_lookback(
152		const arma::mat& conn_counts,
153		const std::vector<std::vector<int>>& unit_row_indices,
154		const std::vector<int>& lookback_sizes,
155		bool aggregate_mean = false,
156		bool weighting_linear = false,
157		const std::vector<int>& segment_ids = {}
158		) {
159	!	arma::mat out = conn_counts;
160	!	size_t n_units = unit_row_indices.size();
161	!	bool has_segments = (!segment_ids.empty() && segment_ids.size() >= conn_counts.n_rows);
162
163	!	for (size_t u = 0; u < n_units; ++u) {
164	!	const auto& rows = unit_row_indices[u];
165	!	if (rows.empty()) continue;
166
167	!	int k = (u < lookback_sizes.size()) ? lookback_sizes[u] : 20;
168	!	int n_u = static_cast<int>(rows.size());
169
170	!	arma::mat block(n_u, conn_counts.n_cols);
171	!	std::vector<int> block_segs(n_u, 0);
172
173	!	for (int i = 0; i < n_u; ++i) {
174	!	int r = rows[i];
175	!	block.row(i) = conn_counts.row(r);
176	!	if (has_segments) block_segs[i] = segment_ids[r];
177		}
178
179		arma::mat smoothed = lookback_block(block, k, aggregate_mean, weighting_linear,
180	!	has_segments ? block_segs : std::vector<int>{});
181
182	!	for (int i = 0; i < n_u; ++i) {
183	!	out.row(rows[i]) = smoothed.row(i);
184		}
185		}
186
187	!	return out;
188		}
189
190		/** @brief Apply EMA door smoothing to an entire connection matrix across multiple units.
191		*/
192	!	inline arma::mat apply_door_ema(
193		const arma::mat& conn_counts,
194		const std::vector<std::vector<int>>& unit_row_indices,
195		double alpha = 0.1,
196		const std::vector<int>& segment_ids = {}
197		) {
198	!	arma::mat out = conn_counts;
199	!	size_t n_units = unit_row_indices.size();
200	!	bool has_segments = (!segment_ids.empty() && segment_ids.size() >= conn_counts.n_rows);
201
202	!	for (size_t u = 0; u < n_units; ++u) {
203	!	const auto& rows = unit_row_indices[u];
204	!	if (rows.empty()) continue;
205
206	!	int n_u = static_cast<int>(rows.size());
207	!	arma::mat block(n_u, conn_counts.n_cols);
208	!	std::vector<int> block_segs(n_u, 0);
209
210	!	for (int i = 0; i < n_u; ++i) {
211	!	int r = rows[i];
212	!	block.row(i) = conn_counts.row(r);
213	!	if (has_segments) block_segs[i] = segment_ids[r];
214		}
215
216	!	arma::mat smoothed = ema_block(block, alpha, has_segments ? block_segs : std::vector<int>{});
217
218	!	for (int i = 0; i < n_u; ++i) {
219	!	out.row(rows[i]) = smoothed.row(i);
220		}
221		}
222
223	!	return out;
224		}
225
226		} // namespace qe
227
228		#endif // LIBQE_DOOR_HPP

1		/** @file linalg_fallback.hpp
2		* @brief Thin LAPACK-free implementations of SVD, QR, and symmetric
3		* eigendecomposition for use in WASM / no-LAPACK builds.
4		*
5		* When @c ARMA_USE_LAPACK is defined these wrappers delegate directly to the
6		* standard Armadillo routines (@c arma::svd, @c arma::qr, @c arma::eig_sym),
7		* so there is zero overhead in normal builds. When LAPACK is not available
8		* (e.g. Emscripten WASM builds) the wrappers fall back to pure C++
9		* implementations that work for the small matrices typical in ENA
10		* (p ≤ ~100 connection dimensions).
11		*
12		* All routines live in namespace @c qe::linalg and operate on @c arma::mat
13		* (double precision).
14		*
15		* Fallback algorithms:
16		* - SVD — one-sided Jacobi iterations on \f$A^T\f$ (Demmel & Veselić 1992)
17		* - QR — modified Gram-Schmidt
18		* - eig_sym — cyclic Jacobi (classical Jacobi sweep)
19		*/
20
21		#ifndef LIBQE_LINALG_FALLBACK_HPP
22		#define LIBQE_LINALG_FALLBACK_HPP
23
24		#include <armadillo>
25		#include <cmath>
26		#include <limits>
27
28		namespace qe {
29		namespace linalg {
30
31		/// @cond INTERNAL
32		// ── helpers ───────────────────────────────────────────────────────────────────
33
34		namespace detail {
35
36		/// @brief Modified Gram-Schmidt orthonormalization applied in-place.
37		///
38		/// Columns of @p Q are successively orthonormalized against all previously
39		/// processed columns. A column whose norm falls below 1e-14 after projection
40		/// is left unchanged (not zeroed) to avoid introducing NaNs.
41		///
42		/// @param Q Matrix whose columns are orthonormalized in-place (modified).
43		inline void mgs(arma::mat& Q) {
44		const arma::uword p = Q.n_cols;
45		for (arma::uword j = 0; j < p; ++j) {
46		for (arma::uword i = 0; i < j; ++i)
47		Q.col(j) -= arma::dot(Q.col(i), Q.col(j)) * Q.col(i);
48		double n = arma::norm(Q.col(j));
49		if (n > 1e-14) Q.col(j) /= n;
50		}
51		}
52
53		/// @brief One cyclic Jacobi sweep on a p×p symmetric matrix.
54		///
55		/// Iterates over all super-diagonal index pairs @c (q,r) and applies a
56		/// two-sided Jacobi rotation that annihilates @c A(q,r). The same rotation
57		/// is accumulated in @p V so that, starting with @c V = I, the full sequence
58		/// of sweeps yields @c V = eigenvector matrix.
59		///
60		/// @param A Symmetric p×p matrix, updated in-place toward diagonal form.
61		/// @param V Rotation accumulator (p×p); initialize to identity to recover
62		/// eigenvectors.
63		///
64		/// @note Algorithm follows the classical cyclic-by-rows Jacobi scheme.
65		/// Off-diagonal entries no larger than 1e-14 × (\|a_qq\| + \|a_rr\|) are
66		/// skipped for numerical stability. The `<=` (not `<`) guard correctly
67		/// handles the all-zero case (a_qq = a_rr = a_qr = 0) that would
68		/// otherwise produce NaN via 0/0 in the theta computation.
69		inline void jacobi_sweep(arma::mat& A, arma::mat& V) {
70		const arma::uword p = A.n_rows;
71		for (arma::uword q = 0; q < p - 1; ++q) {
72		for (arma::uword r = q + 1; r < p; ++r) {
73		double aqq = A(q, q), arr = A(r, r), aqr = A(q, r);
74		if (std::abs(aqr) <= 1e-14 * (std::abs(aqq) + std::abs(arr)))
75		continue;
76		double theta = (arr - aqq) / (2.0 * aqr);
77		double t = (theta >= 0.0)
78		? 1.0 / (theta + std::sqrt(1.0 + theta * theta))
79		: 1.0 / (theta - std::sqrt(1.0 + theta * theta));
80		double c = 1.0 / std::sqrt(1.0 + t * t); ///< cosine of Jacobi rotation angle
81		double s = t * c; ///< sine of Jacobi rotation angle
82		// Apply Jacobi rotation to A (symmetric update)
83		double new_qq = aqq - t * aqr;
84		double new_rr = arr + t * aqr;
85		A(q, q) = new_qq; A(r, r) = new_rr; A(q, r) = A(r, q) = 0.0;
86		for (arma::uword k = 0; k < p; ++k) {
87		if (k == q \|\| k == r) continue;
88		double akq = A(k, q), akr = A(k, r);
89		A(k, q) = A(q, k) = c * akq - s * akr;
90		A(k, r) = A(r, k) = s * akq + c * akr;
91		}
92		// Accumulate rotation in V
93		for (arma::uword k = 0; k < p; ++k) {
94		double vkq = V(k, q), vkr = V(k, r);
95		V(k, q) = c * vkq - s * vkr;
96		V(k, r) = s * vkq + c * vkr;
97		}
98		}
99		}
100		}
101
102		} // namespace detail
103
104		/// @endcond
105
106		/// @name Eigendecomposition
107		/// @{
108
109		// ── eig_sym ───────────────────────────────────────────────────────────────────
110
111		/// @brief Eigendecomposition of a real symmetric matrix.
112		///
113		/// Computes eigenvalues and eigenvectors of the symmetric p×p matrix @p S and
114		/// returns them sorted in ascending eigenvalue order.
115		///
116		/// @param[out] eigval p×1 vector of eigenvalues, ascending.
117		/// @param[out] eigvec p×p matrix whose columns are the corresponding
118		/// orthonormal eigenvectors.
119		/// @param[in] S Real symmetric p×p input matrix (not modified).
120		///
121		/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::eig_sym
122		/// directly. Without LAPACK the fallback runs up to @c 50*p cyclic
123		/// Jacobi sweeps; convergence is guaranteed for well-conditioned
124		/// matrices but may be slow for p > 200.
125		///
126		/// @note Fallback algorithm: classical cyclic-by-rows Jacobi iteration.
127		/// Convergence is declared when the Frobenius norm of the strictly
128		/// off-diagonal part drops below 1e-13.
129	1x	inline void eig_sym(arma::vec& eigval, arma::mat& eigvec, const arma::mat& S) {
130		#if defined(ARMA_USE_LAPACK)
131	1x	arma::eig_sym(eigval, eigvec, S);
132		#else
133		const arma::uword p = S.n_rows;
134		arma::mat A = S; ///< working copy of S
135		eigvec.eye(p, p); ///< accumulate rotations; starts as identity
136
137		for (int sweep = 0; sweep < 50 * static_cast<int>(p); ++sweep) {
138		// Check off-diagonal convergence
139		double off = 0.0;
140		for (arma::uword q = 0; q < p - 1; ++q)
141		for (arma::uword r = q + 1; r < p; ++r)
142		off += A(q, r) * A(q, r);
143		if (std::sqrt(off) < 1e-13) break;
144		detail::jacobi_sweep(A, eigvec);
145		}
146		eigval = A.diag();
147
148		// Sort ascending
149		arma::uvec idx = arma::sort_index(eigval);
150		eigval = eigval(idx);
151		eigvec = eigvec.cols(idx);
152		#endif
153		}
154
155		/// @}
156
157		/// @name QR Factorization
158		/// @{
159
160		// ── qr_full ───────────────────────────────────────────────────────────────────
161
162		/// @brief Full (square) QR factorization of an n×p matrix.
163		///
164		/// Produces a square n×n orthogonal factor @p Q and an n×p upper-trapezoidal
165		/// factor @p R such that @c A = Q*R. This matches the convention of
166		/// @c arma::qr(), which returns the full (non-economy) @p Q.
167		///
168		/// @param[out] Q n×n orthogonal matrix.
169		/// @param[out] R n×p upper-trapezoidal matrix.
170		/// @param[in] A n×p input matrix (not modified).
171		///
172		/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::qr directly.
173		/// Without LAPACK the fallback constructs the full ONB by seeding
174		/// columns beyond @c A.n_cols with standard basis vectors and then
175		/// applying modified Gram-Schmidt; numerical quality degrades for
176		/// nearly rank-deficient inputs.
177		///
178		/// @note The caller (@c orthogonal_svd) depends on @p Q being n×n so that
179		/// @c Q.cols(k, n-1) spans the null space of @c A^T.
180		inline void qr_full(arma::mat& Q, arma::mat& R, const arma::mat& A) {
181		#if defined(ARMA_USE_LAPACK)
182		arma::qr(Q, R, A);
183		#else
184		// Build a full n×n orthogonal matrix.
185		// Step 1: thin QR via modified Gram-Schmidt (gives first k=A.n_cols columns).
186		// Step 2: extend to a full ONB using random vectors + re-orthogonalisation.
187		const arma::uword n = A.n_rows, p = A.n_cols;
188
189		// Start with all n columns of A padded with random vectors if n > p.
190		arma::mat Qfull(n, n, arma::fill::zeros);
191		// Copy A's columns first
192		for (arma::uword j = 0; j < p && j < n; ++j)
193		Qfull.col(j) = A.col(j);
194		// Fill remaining columns with standard basis vectors (or random fallback)
195		for (arma::uword j = p; j < n; ++j) {
196		Qfull(j, j) = 1.0; ///< standard basis vector e_j used as seed
197		}
198		// Full MGS orthonormalization
199		detail::mgs(Qfull);
200		Q = Qfull;
201
202		// R is n×p (upper trapezoidal); only first p rows are non-zero
203		R = Q.t() * A;
204		// Zero the strictly lower-triangular part of R[0:p, 0:p]
205		for (arma::uword j = 0; j < p; ++j)
206		for (arma::uword i = j + 1; i < n; ++i)
207		R(i, j) = 0.0;
208		#endif
209		}
210
211		/// @}
212
213		/// @name Singular Value Decomposition
214		/// @{
215
216		// ── svd ───────────────────────────────────────────────────────────────────────
217
218		/// @brief Full SVD: @c A ≈ U * diag(s) * V', singular values descending.
219		///
220		/// Decomposes the n×p matrix @p A into left singular vectors @p U, singular
221		/// values @p s, and right singular vectors @p V.
222		///
223		/// @param[out] U n×min(n,p) matrix of left singular vectors.
224		/// @param[out] s min(n,p)×1 vector of singular values, descending.
225		/// @param[out] V p×p matrix of right singular vectors (all p columns).
226		/// @param[in] A n×p input matrix (not modified).
227		///
228		/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::svd directly.
229		/// Without LAPACK the fallback computes @c A^T*A, calls @c eig_sym,
230		/// and derives @p U via @c U[:,j] = A*V[:,j]/s_j. Columns with
231		/// @c s_j < 1e-14 are zeroed. Accuracy is governed by the
232		/// conditioning of @c A^T*A (squares the condition number of @p A).
233		///
234		/// @note Fallback algorithm: eigendecomposition of \f$A^T A\f$ via cyclic
235		/// Jacobi (Demmel & Veselić 1992). Suitable for the small, moderately
236		/// conditioned matrices typical in ENA (p ≤ ~100).
237		inline void svd(arma::mat& U, arma::vec& s, arma::mat& V, const arma::mat& A) {
238		#if defined(ARMA_USE_LAPACK)
239		arma::svd(U, s, V, A);
240		#else
241		// One-sided Jacobi on A^T:
242		// We work with B = A^T (p×n). The right singular vectors of A are the
243		// left singular vectors of B (= columns of U_B after Jacobi on B^T B).
244		// We use the eigendecomposition of A^T * A to get V and s.
245		const arma::uword n = A.n_rows, p = A.n_cols;
246		const arma::uword k = std::min(n, p);
247
248		arma::mat AtA = A.t() * A; ///< p×p symmetric Gram matrix
249		arma::vec eigval;
250		arma::mat eigvec;
251		eig_sym(eigval, eigvec, AtA);
252
253		// Eigenvalues from eig_sym are ascending; we want descending.
254		arma::uvec idx = arma::sort_index(eigval, "descend");
255		eigval = eigval(idx);
256		eigvec = eigvec.cols(idx);
257
258		// Singular values (non-negative square roots)
259		s.set_size(k);
260		for (arma::uword j = 0; j < k; ++j)
261		s(j) = std::sqrt(std::max(0.0, eigval(j)));
262
263		// V = eigvec (right singular vectors of A are eigvectors of A^T*A)
264		V = eigvec;
265
266		// U = A * V / s (left singular vectors, only first k columns needed)
267		U.set_size(n, k);
268		for (arma::uword j = 0; j < k; ++j) {
269		if (s(j) > 1e-14)
270		U.col(j) = A * V.col(j) / s(j);
271		else
272		U.col(j).zeros();
273		}
274		#endif
275		}
276
277		// ── svd_right ─────────────────────────────────────────────────────────────────
278
279		/// @brief Singular values and right singular vectors only (left vectors skipped).
280		///
281		/// A lighter-weight variant of @c svd() for callers that only need @p s and
282		/// @p V. Under LAPACK the full @c arma::svd is still called internally; the
283		/// savings come in WASM builds where computing @p U is avoided.
284		///
285		/// @param[out] s min(n,p)×1 vector of singular values, descending.
286		/// @param[out] V p×p matrix of right singular vectors.
287		/// @param[in] A n×p input matrix (not modified).
288		///
289		/// @warning With @c ARMA_USE_LAPACK defined this still invokes the full
290		/// @c arma::svd; the left-vector output is discarded.
291		inline void svd_right(arma::vec& s, arma::mat& V, const arma::mat& A) {
292		#if defined(ARMA_USE_LAPACK)
293		arma::mat U;
294		arma::svd(U, s, V, A);
295		#else
296		arma::mat U_unused;
297		svd(U_unused, s, V, A);
298		#endif
299		}
300
301		/// @}
302
303		/// @name Linear System Solve
304		/// @{
305
306		// ── solve_spd ─────────────────────────────────────────────────────────────────
307
308		/// @brief Solve @c A*x = b for symmetric positive (semi-)definite @p A.
309		///
310		/// Returns the n×k solution matrix @p x given the n×n system matrix @p A and
311		/// the n×k right-hand side @p b. Near-zero eigenvalues are handled gracefully,
312		/// making this a pseudo-inverse solve safe for rank-deficient Gram matrices.
313		///
314		/// @param[in] A Real symmetric positive (semi-)definite n×n matrix.
315		/// @param[in] b n×k right-hand side matrix.
316		/// @param[in] tol Relative threshold below which eigenvalues are treated as
317		/// zero. Eigenvalues with @c \|λ\| < tol * max(\|λ\|) are
318		/// inverted as zero (i.e. those directions are projected out).
319		/// Default: 1e-12.
320		/// @returns n×k solution matrix @p x satisfying @c A*x ≈ b.
321		///
322		/// @warning With @c ARMA_USE_LAPACK defined this delegates to
323		/// @c arma::solve(A, b, arma::solve_opts::fast), which uses a
324		/// Cholesky factorization and does not handle rank deficiency.
325		/// Without LAPACK the eigendecomposition-based pseudo-inverse is used,
326		/// which is safe for rank-deficient inputs but slower.
327		///
328		/// @note Without LAPACK: @c A = V * diag(λ) * V^T is computed via @c eig_sym;
329		/// the solution is @c x = V * diag(1/λ) * V^T * b with small λ skipped.
330		inline arma::mat solve_spd(const arma::mat& A, const arma::mat& b,
331		double tol = 1e-12) {
332		#if defined(ARMA_USE_LAPACK)
333		return arma::solve(A, b, arma::solve_opts::fast);
334		#else
335		arma::vec eigval;
336		arma::mat eigvec;
337		eig_sym(eigval, eigvec, A); ///< A = eigvec * diag(eigval) * eigvec^T
338
339		const double thresh = tol * std::abs(eigval.max());
340		arma::vec inv_eigval(eigval.n_elem, arma::fill::zeros);
341		for (arma::uword i = 0; i < eigval.n_elem; ++i)
342		if (std::abs(eigval(i)) > thresh)
343		inv_eigval(i) = 1.0 / eigval(i);
344
345		// x = eigvec * diag(inv_eigval) * eigvec^T * b
346		return eigvec * arma::diagmat(inv_eigval) * eigvec.t() * b;
347		#endif
348		}
349
350		// ── inv_sym ───────────────────────────────────────────────────────────────────
351
352		/// @brief Inverse of a symmetric positive (semi-)definite matrix, falling back
353		/// to the pseudo-inverse when it is singular.
354		///
355		/// @param[in] A Real symmetric positive (semi-)definite n×n matrix.
356		/// @param[in] tol Without LAPACK: relative eigenvalue threshold, as in
357		/// @c solve_spd. Default: 1e-12.
358		/// @returns n×n matrix: @c A^-1, or the pseudo-inverse when @p A is singular.
359		///
360		/// @warning With @c ARMA_USE_LAPACK defined this calls @c arma::inv and, when
361		/// that reports a singular matrix, @c arma::pinv. Without LAPACK it
362		/// is the eigendecomposition-based pseudo-inverse (@c solve_spd with
363		/// the identity), which equals the inverse for a non-singular @p A.
364	5x	inline arma::mat inv_sym(const arma::mat& A, double tol = 1e-12) {
365		#if defined(ARMA_USE_LAPACK)
366	5x	arma::mat out;
367	5x	if (!arma::inv(out, A))
368	!	out = arma::pinv(A);
369	5x	return out;
370		#else
371		return solve_spd(A, arma::eye(A.n_rows, A.n_rows), tol);
372		#endif
373		}
374
375		/// @}
376
377		} // namespace linalg
378		} // namespace qe
379
380		#endif // LIBQE_LINALG_FALLBACK_HPP

1		/**
2		* @file normalization.hpp
3		* @brief Row-wise normalization utilities for ENA connection networks.
4		*/
5		#ifndef LIBQE_NORMALIZATION_HPP
6		#define LIBQE_NORMALIZATION_HPP
7
8		#include <armadillo>
9		#include <cmath>
10
11		namespace qe {
12
13		/**
14		* @brief Normalize each row to unit L2 norm (project onto the unit hypersphere).
15		*
16		* Rows whose L2 norm is zero are left unchanged.
17		*
18		* @param m Input matrix (n_units × n_connections), row-major.
19		* @returns Matrix of the same shape with each row divided by its L2 norm.
20		* @note Equivalent to @c fun_sphere_norm() in rENA/ena.cpp.
21		*/
22	4x	inline arma::mat normalize_networks(arma::mat m) {
23	4x	arma::mat out(m.n_rows, m.n_cols, arma::fill::zeros);
24	14x	for (arma::uword r = 0; r < m.n_rows; r++) {
25	10x	double len = arma::norm(m.row(r), 2);
26	18x	if (std::isfinite(len) && len > 0.0) out.row(r) = m.row(r) / len;
27		}
28	4x	return out;
29		}
30
31		/**
32		* @brief Scale every entry by the reciprocal of the largest row L2 norm.
33		*
34		* Preserves relative magnitudes across rows; does not normalise each row
35		* individually. If the largest row norm is zero, the matrix is returned
36		* unchanged.
37		*
38		* @param m Input matrix (n_units × n_connections), modified in place.
39		* @returns Scaled matrix (same shape as @p m).
40		* @note Equivalent to @c fun_skip_sphere_norm() in rENA/ena.cpp.
41		*/
42	3x	inline arma::mat scale_networks(arma::mat m) {
43	3x	double largest = 0.0;
44	10x	for (arma::uword r = 0; r < m.n_rows; r++) {
45	7x	double len = arma::norm(m.row(r), 2);
46	7x	if (std::isfinite(len)) largest = std::max(largest, len);
47		}
48	3x	if (largest > 0.0) m /= largest;
49	6x	return m;
50		}
51
52		} // namespace qe
53
54		#endif // LIBQE_NORMALIZATION_HPP

1		/** @file stability.hpp
2		* @brief Stability metrics and distance-distance matrix correlation for parameter tuning.
3		*/
4		#ifndef LIBQE_STABILITY_HPP
5		#define LIBQE_STABILITY_HPP
6
7		#include <armadillo>
8		#include <cmath>
9		#include <algorithm>
10
11		namespace qe {
12
13		/** @brief Compute Pearson correlation between pairwise Euclidean distance structures of two matrices.
14		*
15		* Evaluates dist-dist correlation: r(dist(X), dist(Y)), measuring whether relative
16		* neighborhoods and distances are preserved across parameter changes (e.g. window vs lookback).
17		*
18		* @param[in] X Numeric matrix (n_rows × n_cols).
19		* @param[in] Y Numeric matrix (n_rows × n_cols).
20		*
21		* @returns Pearson correlation coefficient in [-1, 1].
22		*/
23	1x	inline double dist_dist_correlation(const arma::mat& X, const arma::mat& Y) {
24	1x	size_t n = X.n_rows;
25		if (n != Y.n_rows \|\| n < 3) return 0.0;
26
27	1x	size_t n_pairs = n * (n - 1) / 2;
28	1x	arma::vec dx(n_pairs);
29	1x	arma::vec dy(n_pairs);
30
31	1x	size_t idx = 0;
32	11x	for (size_t i = 0; i < n; ++i) {
33	55x	for (size_t j = i + 1; j < n; ++j) {
34	180x	arma::rowvec diff_x = X.row(i) - X.row(j);
35	180x	arma::rowvec diff_y = Y.row(i) - Y.row(j);
36	90x	dx[idx] = std::sqrt(arma::dot(diff_x, diff_x));
37	45x	dy[idx] = std::sqrt(arma::dot(diff_y, diff_y));
38	45x	idx++;
39		}
40		}
41
42	1x	double mean_x = arma::mean(dx);
43	1x	double mean_y = arma::mean(dy);
44
45	1x	arma::vec zx = dx - mean_x;
46	1x	arma::vec zy = dy - mean_y;
47
48	1x	double sxx = arma::dot(zx, zx);
49	1x	double syy = arma::dot(zy, zy);
50	1x	double sxy = arma::dot(zx, zy);
51
52		if (sxx <= 1e-12 \|\| syy <= 1e-12) return 0.0;
53	1x	return sxy / std::sqrt(sxx * syy);
54		}
55
56		} // namespace qe
57
58		#endif // LIBQE_STABILITY_HPP

1		/** @file trajectory_distance.hpp
2		* @brief Integrated Euclidean distance between trajectory curves and lagged distance metrics.
3		*/
4		#ifndef LIBQE_TRAJECTORY_DISTANCE_HPP
5		#define LIBQE_TRAJECTORY_DISTANCE_HPP
6
7		#include <armadillo>
8		#include <vector>
9		#include <cmath>
10		#include <algorithm>
11		#include "trajectory.hpp"
12		#include "linalg_fallback.hpp"
13
14		namespace qe {
15
16		struct GaussLegendreRule {
17		arma::vec nodes;
18		arma::vec weights;
19		};
20
21		/** @brief Obtain exact 32-point Gauss-Legendre quadrature nodes and weights on [-1, 1] via Golub-Welsch. */
22	2x	inline const GaussLegendreRule& get_gl_rule_32() {
23	1x	static const GaussLegendreRule rule = []() {
24	1x	int N = 32;
25	1x	arma::mat J(N, N, arma::fill::zeros);
26	32x	for (int k = 1; k < N; ++k) {
27	31x	double b = static_cast<double>(k) / std::sqrt(4.0 * k * k - 1.0);
28	31x	J(k, k - 1) = b;
29	62x	J(k - 1, k) = b;
30		}
31	1x	arma::vec evals;
32	1x	arma::mat evecs;
33	1x	qe::linalg::eig_sym(evals, evecs, J); // LAPACK-free in WASM builds
34	1x	arma::vec w = 2.0 * arma::square(evecs.row(0).t());
35		return GaussLegendreRule{evals, w};
36		}();
37	2x	return rule;
38		}
39
40		/** @brief Compute integrated Euclidean distance between two 2D parametric polynomial curves.
41		*
42		* Computes integral_{t_start}^{t_end} \|\|F_a(t) - F_b(t)\|\| dt
43		* using 32-point Gauss-Legendre quadrature.
44		*
45		* @param[in] coeffs_ax Coefficients of x_a(t).
46		* @param[in] coeffs_ay Coefficients of y_a(t).
47		* @param[in] coeffs_bx Coefficients of x_b(t).
48		* @param[in] coeffs_by Coefficients of y_b(t).
49		* @param[in] t_start Start of integration interval (default 0.0).
50		* @param[in] t_end End of integration interval (default 1.0).
51		*
52		* @returns Integrated distance scalar.
53		*/
54	2x	inline double integrated_curve_distance(
55		const arma::vec& coeffs_ax,
56		const arma::vec& coeffs_ay,
57		const arma::vec& coeffs_bx,
58		const arma::vec& coeffs_by,
59		double t_start = 0.0,
60		double t_end = 1.0
61		) {
62		if (t_end <= t_start) return 0.0;
63
64	2x	const auto& rule = get_gl_rule_32();
65	2x	double mid = 0.5 * (t_start + t_end);
66	2x	double half_len = 0.5 * (t_end - t_start);
67	2x	double total = 0.0;
68
69	66x	for (size_t i = 0; i < rule.nodes.n_elem; ++i) {
70	64x	double t_eval = mid + half_len * rule.nodes[i];
71	64x	double xa = eval_poly_scalar(coeffs_ax, t_eval);
72	64x	double ya = eval_poly_scalar(coeffs_ay, t_eval);
73	64x	double xb = eval_poly_scalar(coeffs_bx, t_eval);
74	64x	double yb = eval_poly_scalar(coeffs_by, t_eval);
75	64x	double dx = xa - xb;
76	64x	double dy = ya - yb;
77	64x	double dist = std::sqrt(dx * dx + dy * dy);
78
79	128x	total += rule.weights[i] * dist;
80		}
81
82	2x	return half_len * total;
83		}
84
85		/** @brief Compute lagged curve distance: 1/(1-lag) * integral_{lag}^1 \|\|F_follower(t) - F_leader(t - lag)\|\| dt. */
86	!	inline double lagged_curve_distance(
87		const arma::vec& coeffs_fol_x,
88		const arma::vec& coeffs_fol_y,
89		const arma::vec& coeffs_ldr_x,
90		const arma::vec& coeffs_ldr_y,
91		double lag = 0.0
92		) {
93	!	if (lag < 0.0 \|\| lag >= 0.999) return 0.0;
94	!	double t_start = lag;
95	!	double t_end = 1.0;
96	!	const auto& rule = get_gl_rule_32();
97	!	double mid = 0.5 * (t_start + t_end);
98	!	double half_len = 0.5 * (t_end - t_start);
99	!	double total = 0.0;
100
101	!	for (size_t i = 0; i < rule.nodes.n_elem; ++i) {
102	!	double t_fol = mid + half_len * rule.nodes[i];
103	!	double t_ldr = t_fol - lag;
104
105	!	double x_fol = eval_poly_scalar(coeffs_fol_x, t_fol);
106	!	double y_fol = eval_poly_scalar(coeffs_fol_y, t_fol);
107	!	double x_ldr = eval_poly_scalar(coeffs_ldr_x, t_ldr);
108	!	double y_ldr = eval_poly_scalar(coeffs_ldr_y, t_ldr);
109
110	!	double dx = x_fol - x_ldr;
111	!	double dy = y_fol - y_ldr;
112	!	double dist = std::sqrt(dx * dx + dy * dy);
113
114	!	total += rule.weights[i] * dist;
115		}
116
117	!	double integral = half_len * total;
118	!	double interval = 1.0 - lag;
119	!	return (interval > 1e-6) ? (integral / interval) : 0.0;
120		}
121
122		/** @brief Compute pairwise integrated distance matrix across K polynomial trajectory curves. */
123	!	inline arma::mat pairwise_trajectory_distance_matrix(
124		const std::vector<arma::vec>& all_coeffs_x,
125		const std::vector<arma::vec>& all_coeffs_y
126		) {
127	!	size_t K = all_coeffs_x.size();
128	!	arma::mat D(K, K, arma::fill::zeros);
129
130	!	for (size_t i = 0; i < K; ++i) {
131	!	for (size_t j = i + 1; j < K; ++j) {
132	!	double dist = integrated_curve_distance(
133	!	all_coeffs_x[i], all_coeffs_y[i],
134	!	all_coeffs_x[j], all_coeffs_y[j]
135		);
136	!	D(i, j) = dist;
137	!	D(j, i) = dist;
138		}
139		}
140
141	!	return D;
142		}
143
144		} // namespace qe
145
146		#endif // LIBQE_TRAJECTORY_DISTANCE_HPP

1		/** @file trajectory_following.hpp
2		* @brief Who-follows-whom leader-follower turn lag analysis along shared timeline.
3		*/
4		#ifndef LIBQE_TRAJECTORY_FOLLOWING_HPP
5		#define LIBQE_TRAJECTORY_FOLLOWING_HPP
6
7		#include <armadillo>
8		#include <vector>
9		#include <cmath>
10		#include <algorithm>
11		#include <limits>
12
13		namespace qe {
14
15		/** @brief Result of signed turn-lag sweep across delta in [-max_lag, max_lag]. */
16		struct SignedTurnLagResult {
17		int best_lag = 0; ///< delta* that minimizes mean distance
18		double min_mean_distance = 0.0;
19		arma::ivec lags; ///< Range of evaluated lags
20		arma::vec mean_distances; ///< Mean distance for each lag
21		arma::ivec valid_counts; ///< Number of paired turns for each lag
22		};
23
24		/** @brief Compute mean distance between agent A at turn T and agent B at turn T - delta.
25		*
26		* @param[in] pts_a n_a × 2 coordinates for agent A.
27		* @param[in] pts_b n_b × 2 coordinates for agent B.
28		* @param[in] times_a Length n_a vector of calendar/turn indices for A.
29		* @param[in] times_b Length n_b vector of calendar/turn indices for B.
30		* @param[in] delta Signed lag in turns (delta > 0: A compares to B's past => A follows B;
31		* delta < 0: B compares to A's past => B follows A; delta = 0: sync).
32		*
33		* @returns Pair of (mean_distance, valid_pairs_count).
34		*/
35	9x	inline std::pair<double, int> signed_turn_lag_distance(
36		const arma::mat& pts_a,
37		const arma::mat& pts_b,
38		const arma::vec& times_a,
39		const arma::vec& times_b,
40		int delta
41		) {
42	9x	int n_a = static_cast<int>(pts_a.n_rows);
43	9x	int n_b = static_cast<int>(pts_b.n_rows);
44		if (n_a == 0 \|\| n_b == 0) return {0.0, 0};
45
46	9x	double sum_dist = 0.0;
47	9x	int count = 0;
48
49	54x	for (int i = 0; i < n_a; ++i) {
50	45x	double t_a = times_a[i];
51	45x	double target_tb = t_a - static_cast<double>(delta);
52
53		// Find the latest turn for agent B that occurred at or before target_tb
54	45x	int best_b_idx = -1;
55	216x	for (int j = 0; j < n_b; ++j) {
56	382x	if (times_b[j] <= target_tb) {
57	171x	best_b_idx = j;
58		} else {
59	20x	break;
60		}
61		}
62
63	45x	if (best_b_idx >= 0) {
64	84x	double dx = pts_a(i, 0) - pts_b(best_b_idx, 0);
65	84x	double dy = pts_a(i, 1) - pts_b(best_b_idx, 1);
66	42x	sum_dist += std::sqrt(dx * dx + dy * dy);
67	42x	count++;
68		}
69		}
70
71		double mean_d = (count > 0) ? (sum_dist / count) : 0.0;
72	9x	return {mean_d, count};
73		}
74
75		/** @brief Sweep signed turn lags in [-max_lag, max_lag] to determine leader-follower dynamics. */
76	1x	inline SignedTurnLagResult best_signed_turn_lag(
77		const arma::mat& pts_a,
78		const arma::mat& pts_b,
79		const arma::vec& times_a,
80		const arma::vec& times_b,
81		int max_lag = 15
82		) {
83	1x	int k_max = std::max(1, max_lag);
84	1x	int n_lags = 2 * k_max + 1;
85
86	1x	SignedTurnLagResult res;
87	1x	res.lags.set_size(n_lags);
88	1x	res.mean_distances.set_size(n_lags);
89	1x	res.valid_counts.set_size(n_lags);
90
91	1x	double min_dist = std::numeric_limits<double>::infinity();
92	1x	int best_delta = 0;
93
94	10x	for (int idx = 0; idx < n_lags; ++idx) {
95	9x	int delta = -k_max + idx;
96	9x	res.lags[idx] = delta;
97
98	9x	auto [mean_d, count] = signed_turn_lag_distance(pts_a, pts_b, times_a, times_b, delta);
99	9x	res.mean_distances[idx] = mean_d;
100	9x	res.valid_counts[idx] = count;
101
102	9x	if (count > 0 && mean_d < min_dist) {
103	5x	min_dist = mean_d;
104	5x	best_delta = delta;
105		}
106		}
107
108	1x	res.best_lag = best_delta;
109		res.min_mean_distance = (min_dist == std::numeric_limits<double>::infinity()) ? 0.0 : min_dist;
110	1x	return res;
111		}
112
113		} // namespace qe
114
115		#endif // LIBQE_TRAJECTORY_FOLLOWING_HPP

1		/** @file trajectory.hpp
2		* @brief Parametric trajectory curve fitting, LOOCV-2D degree selection,
3		* and differential geometry (velocity, acceleration, curvature, turns).
4		*/
5		#ifndef LIBQE_TRAJECTORY_HPP
6		#define LIBQE_TRAJECTORY_HPP
7
8		#include <armadillo>
9		#include <vector>
10		#include <string>
11		#include <cmath>
12		#include <algorithm>
13		#include <limits>
14
15		#include "linalg_fallback.hpp"
16
17		namespace qe {
18
19		/** @brief Polynomial curve fit result in 2D epistemic space. */
20		struct PolyCurveFit {
21		int degree = 1;
22		arma::vec coeffs_x; ///< Raw coefficients [c0, ..., cd] for x(t) = sum(ck * t^k)
23		arma::vec coeffs_y; ///< Raw coefficients [c0, ..., cd] for y(t) = sum(ck * t^k)
24		std::string basis = "orthogonal"; ///< Basis used during fitting: "orthogonal" (R stats::poly-like) or "raw".
25		arma::vec basis_alpha; ///< Orthogonal-basis recurrence alpha values (empty for raw basis).
26		arma::vec basis_norm2; ///< Orthogonal-basis recurrence squared norms, matching stats::poly coefs$norm2.
27		double cv_error = 0.0; ///< 2D LOOCV error
28		double aic = 0.0; ///< Joint bivariate AIC
29		arma::vec t; ///< Evaluated or input time points
30		arma::mat fitted; ///< n × 2 fitted coordinates
31		};
32
33	5x	inline arma::mat build_raw_poly_design(const arma::vec& t_vec, int degree) {
34	5x	int n = static_cast<int>(t_vec.n_elem);
35	5x	arma::mat X(n, degree + 1);
36	50x	for (int i = 0; i < n; ++i) {
37	45x	double cur_t = t_vec[i];
38	45x	double p = 1.0;
39	183x	for (int k = 0; k <= degree; ++k) {
40	138x	X(i, k) = p;
41	138x	p *= cur_t;
42		}
43		}
44	5x	return X;
45		}
46
47	5x	inline arma::vec poly_multiply_by_linear(const arma::vec& coeffs, double shift) {
48	5x	arma::vec out(coeffs.n_elem + 1, arma::fill::zeros);
49	17x	for (size_t k = 0; k < coeffs.n_elem; ++k) {
50	24x	out[k] += -shift * coeffs[k];
51	36x	out[k + 1] += coeffs[k];
52		}
53	5x	return out;
54		}
55
56	4x	inline arma::mat orthogonal_poly_basis_from_coefs(
57		const arma::vec& t_vec,
58		int degree,
59		const arma::vec& alpha,
60		const arma::vec& norm2
61		) {
62	4x	int n = static_cast<int>(t_vec.n_elem);
63	4x	arma::mat Z(n, degree + 1, arma::fill::ones);
64	4x	if (degree >= 1) {
65	12x	Z.col(1) = t_vec - alpha[0];
66		}
67	9x	for (int i = 2; i <= degree; ++i) {
68	30x	Z.col(i) = (t_vec - alpha[i - 1]) % Z.col(i - 1) -
69	25x	(norm2[i] / norm2[i - 1]) * Z.col(i - 2);
70		}
71	17x	for (int i = 0; i <= degree; ++i) {
72	39x	Z.col(i) /= std::sqrt(norm2[i + 1]);
73		}
74	12x	return Z.cols(1, degree);
75		}
76
77	4x	inline arma::mat orthogonal_basis_raw_coefficients(
78		int degree,
79		const arma::vec& alpha,
80		const arma::vec& norm2
81		) {
82	4x	arma::mat raw(degree + 1, degree + 1, arma::fill::zeros);
83	4x	raw(0, 0) = 1.0;
84	4x	if (degree >= 1) {
85	4x	arma::vec p1(2);
86	12x	p1[0] = -alpha[0];
87	4x	p1[1] = 1.0;
88	12x	raw.submat(0, 1, 1, 1) = p1 / std::sqrt(norm2[2]);
89		}
90	4x	arma::vec prev2(1);
91	4x	prev2[0] = 1.0;
92	4x	arma::vec prev1;
93	4x	if (degree >= 1) {
94	4x	prev1.set_size(2);
95	12x	prev1[0] = -alpha[0];
96	4x	prev1[1] = 1.0;
97		}
98	9x	for (int i = 2; i <= degree; ++i) {
99	10x	arma::vec term = poly_multiply_by_linear(prev1, alpha[i - 1]);
100	5x	arma::vec padded_prev2(term.n_elem, arma::fill::zeros);
101	10x	padded_prev2.subvec(0, prev2.n_elem - 1) = prev2;
102	25x	arma::vec current = term - (norm2[i] / norm2[i - 1]) * padded_prev2;
103	20x	raw.submat(0, i, i, i) = current / std::sqrt(norm2[i + 1]);
104	5x	prev2 = prev1;
105	5x	prev1 = current;
106		}
107	4x	return raw;
108		}
109
110	4x	inline bool modified_gram_schmidt_qr(
111		const arma::mat& X,
112		arma::mat& Q,
113		arma::mat& R
114		) {
115	4x	const int n = static_cast<int>(X.n_rows);
116	4x	const int p = static_cast<int>(X.n_cols);
117	4x	Q.zeros(n, p);
118	4x	R.zeros(p, p);
119
120	17x	for (int k = 0; k < p; ++k) {
121	26x	arma::vec v = X.col(k);
122	29x	for (int j = 0; j < k; ++j) {
123	48x	R(j, k) = arma::dot(Q.col(j), v);
124	64x	v -= R(j, k) * Q.col(j);
125		}
126
127	13x	double norm = std::sqrt(arma::dot(v, v));
128	13x	if (norm <= 1e-12) {
129	!	return false;
130		}
131
132	26x	R(k, k) = norm;
133	26x	Q.col(k) = v / norm;
134		}
135
136	4x	return true;
137		}
138
139	4x	inline arma::mat build_orthogonal_poly_design(
140		const arma::vec& t_vec,
141		int degree,
142		arma::vec& alpha,
143		arma::vec& norm2
144		) {
145	4x	int n = static_cast<int>(t_vec.n_elem);
146	4x	if (degree < 1) {
147	!	return arma::mat(n, 0);
148		}
149	4x	arma::vec unique_t = arma::unique(t_vec);
150	4x	if (static_cast<size_t>(degree) >= unique_t.n_elem) {
151	!	throw std::invalid_argument("orthogonal polynomial degree must be less than number of unique time points.");
152		}
153
154	4x	double xbar = arma::mean(t_vec);
155	4x	arma::vec centered = t_vec - xbar;
156	4x	arma::mat X = build_raw_poly_design(centered, degree);
157
158	4x	arma::mat Q;
159	4x	arma::mat R;
160	4x	bool ok = modified_gram_schmidt_qr(X, Q, R);
161	4x	if (!ok) {
162	!	throw std::invalid_argument("orthogonal polynomial degree must be less than number of unique time points.");
163		}
164
165	8x	arma::mat Z = Q * arma::diagmat(R.diag());
166	4x	arma::vec z_norm2 = arma::sum(arma::square(Z), 0).t();
167
168	4x	alpha.set_size(degree);
169	13x	for (int i = 0; i < degree; ++i) {
170	36x	alpha[i] = arma::dot(centered, arma::square(Z.col(i))) / z_norm2[i] + xbar;
171		}
172
173	4x	norm2.set_size(degree + 2);
174	4x	norm2[0] = 1.0;
175	17x	for (int i = 0; i <= degree; ++i) {
176	39x	norm2[i + 1] = z_norm2[i];
177		}
178
179	4x	return orthogonal_poly_basis_from_coefs(t_vec, degree, alpha, norm2);
180		}
181
182	5x	inline arma::mat build_poly_design(
183		const arma::vec& t_vec,
184		int degree,
185		const std::string& basis,
186		arma::vec& alpha,
187		arma::vec& norm2,
188		arma::mat& basis_raw_coeffs
189		) {
190	5x	int n = static_cast<int>(t_vec.n_elem);
191	5x	arma::mat X(n, degree + 1);
192	5x	X.col(0).ones();
193	5x	if (basis == "orthogonal") {
194	4x	arma::mat Z = build_orthogonal_poly_design(t_vec, degree, alpha, norm2);
195	8x	X.cols(1, degree) = Z;
196	4x	basis_raw_coeffs = orthogonal_basis_raw_coefficients(degree, alpha, norm2);
197	5x	} else if (basis == "raw") {
198	1x	X = build_raw_poly_design(t_vec, degree);
199	1x	alpha.reset();
200	1x	norm2.reset();
201	1x	basis_raw_coeffs.eye(degree + 1, degree + 1);
202		} else {
203	!	throw std::invalid_argument("basis must be 'orthogonal' or 'raw'.");
204		}
205	5x	return X;
206		}
207
208		/** @brief Differential geometry profile along a trajectory curve. */
209		struct TrajectoryDerivatives {
210		arma::vec t;
211		arma::vec vx; ///< dx/dt
212		arma::vec vy; ///< dy/dt
213		arma::vec speed; ///< sqrt(vx^2 + vy^2)
214		arma::vec ax; ///< d^2x/dt^2
215		arma::vec ay; ///< d^2y/dt^2
216		arma::vec heading_rate; ///< omega = vxay - vyax (angular turning velocity)
217		arma::vec curvature; ///< kappa = omega / max(speed^3, eps)
218		};
219
220		/** @brief Evaluate a 1D polynomial sum(c_k * t^k) at a single time t. */
221	296x	inline double eval_poly_scalar(const arma::vec& coeffs, double t) {
222	296x	double val = 0.0;
223	296x	double t_pow = 1.0;
224	956x	for (size_t k = 0; k < coeffs.n_elem; ++k) {
225	660x	val += coeffs[k] * t_pow;
226	660x	t_pow *= t;
227		}
228	296x	return val;
229		}
230
231		/** @brief Evaluate polynomial 1st derivative sum(k * c_k * t^(k-1)) at a single time t. */
232	6x	inline double eval_poly_deriv1_scalar(const arma::vec& coeffs, double t) {
233	6x	double val = 0.0;
234	6x	double t_pow = 1.0;
235	15x	for (size_t k = 1; k < coeffs.n_elem; ++k) {
236	9x	val += static_cast<double>(k) * coeffs[k] * t_pow;
237	9x	t_pow *= t;
238		}
239	6x	return val;
240		}
241
242		/** @brief Evaluate polynomial 2nd derivative sum(k(k-1) c_k * t^(k-2)) at a single time t. */
243	6x	inline double eval_poly_deriv2_scalar(const arma::vec& coeffs, double t) {
244	6x	double val = 0.0;
245	6x	double t_pow = 1.0;
246	9x	for (size_t k = 2; k < coeffs.n_elem; ++k) {
247	3x	val += static_cast<double>(k * (k - 1)) * coeffs[k] * t_pow;
248	3x	t_pow *= t;
249		}
250	6x	return val;
251		}
252
253		/** @brief Evaluate 2D parametric curve (x(t), y(t)) at multiple time points. */
254	2x	inline arma::mat eval_poly_curve(
255		const arma::vec& coeffs_x,
256		const arma::vec& coeffs_y,
257		const arma::vec& t_eval
258		) {
259	2x	size_t n = t_eval.n_elem;
260	2x	arma::mat out(n, 2);
261	22x	for (size_t i = 0; i < n; ++i) {
262	20x	double t = t_eval[i];
263	20x	out(i, 0) = eval_poly_scalar(coeffs_x, t);
264	40x	out(i, 1) = eval_poly_scalar(coeffs_y, t);
265		}
266	2x	return out;
267		}
268
269		/** @brief Fit a 2D polynomial curve with automatic LOOCV-2D or AIC degree selection.
270		*
271		* Given coordinates (x_i, y_i) and times t_i in [0, 1], constructs polynomial models
272		* for degrees d in [1, max_degree]. LOOCV error is computed via the closed-form hat-matrix formula:
273		* CV(d) = sum( (r_{x,i} / (1 - h_{ii}))^2 + (r_{y,i} / (1 - h_{ii}))^2 )
274		*
275		* @param[in] points n × 2 matrix of trajectory coordinates.
276		* @param[in] t Vector of length n of strictly increasing times (default [0, 1] linearly spaced).
277		* @param[in] max_degree Maximum polynomial degree to evaluate (default 3).
278		* @param[in] fixed_degree When >= 1, uses this exact degree without search.
279		* @param[in] criterion "loocv" (default) or "aic".
280		*
281		* @returns PolyCurveFit containing optimal degree, coefficients, errors, and fitted points.
282		*/
283	3x	inline PolyCurveFit fit_poly_loocv2d(
284		const arma::mat& points,
285		const arma::vec& t = arma::vec(),
286		int max_degree = 3,
287		int fixed_degree = 0,
288		const std::string& criterion = "loocv",
289		const std::string& basis = "orthogonal"
290		) {
291	3x	int n = static_cast<int>(points.n_rows);
292	3x	if (n < 2 \|\| points.n_cols < 2) {
293	!	throw std::invalid_argument("fit_poly_loocv2d requires at least 2 points and 2 columns.");
294		}
295
296		// Default normalized time t in [0, 1]
297	3x	arma::vec t_vec;
298	3x	if (t.n_elem == static_cast<size_t>(n)) {
299	3x	t_vec = t;
300		} else {
301	!	t_vec = arma::linspace<arma::vec>(0.0, 1.0, n);
302		}
303
304	3x	arma::vec x = points.col(0);
305	3x	arma::vec y = points.col(1);
306
307	3x	int deg_upper = std::min(max_degree, std::max(1, n - 2));
308	3x	if (fixed_degree >= 1) {
309	2x	deg_upper = std::min(fixed_degree, n - 1);
310		}
311
312	3x	int best_deg = 1;
313	3x	double best_score = std::numeric_limits<double>::infinity();
314	3x	PolyCurveFit best_fit;
315
316	3x	int deg_start = (fixed_degree >= 1) ? fixed_degree : 1;
317
318	8x	for (int d = deg_start; d <= deg_upper; ++d) {
319	5x	int n_params = d + 1;
320		if (n_params > n) break;
321
322	5x	arma::vec basis_alpha;
323	5x	arma::vec basis_norm2;
324	5x	arma::mat basis_raw_coeffs;
325	5x	arma::mat X = build_poly_design(t_vec, d, basis, basis_alpha, basis_norm2, basis_raw_coeffs);
326
327		// Solve least squares: beta = (X^T X)^(-1) X^T y
328	5x	arma::mat XtX = X.t() * X;
329		// Pseudo-inverse when singular; LAPACK-free in WASM builds.
330	5x	arma::mat XtX_inv = qe::linalg::inv_sym(XtX);
331
332	5x	arma::vec beta_x = XtX_inv * (X.t() * x);
333	5x	arma::vec beta_y = XtX_inv * (X.t() * y);
334
335	5x	arma::vec pred_x = X * beta_x;
336	5x	arma::vec pred_y = X * beta_y;
337
338	10x	arma::vec rx = x - pred_x;
339	5x	arma::vec ry = y - pred_y;
340
341		// Hat matrix diagonal: h_ii = x_i^T (X^T X)^(-1) x_i
342	5x	arma::vec h(n);
343	50x	for (int i = 0; i < n; ++i) {
344	90x	arma::rowvec xi = X.row(i);
345	90x	h[i] = arma::as_scalar(xi * XtX_inv * xi.t());
346	90x	h[i] = std::min(0.9999, std::max(0.0, h[i]));
347		}
348
349		// 2D LOOCV error
350	5x	double cv_err = 0.0;
351	50x	for (int i = 0; i < n; ++i) {
352	45x	double denom = 1.0 - h[i];
353	45x	double rxi = rx[i] / denom;
354	45x	double ryi = ry[i] / denom;
355	45x	cv_err += (rxi * rxi + ryi * ryi);
356		}
357
358		// Bivariate AIC
359	5x	double rss = arma::dot(rx, rx) + arma::dot(ry, ry);
360	5x	double aic = static_cast<double>(n) * std::log(std::max(1e-12, rss / (2.0 * n))) + 2.0 * (2.0 * n_params);
361
362		double score = (criterion == "aic") ? aic : cv_err;
363
364	5x	if (score < best_score \|\| d == deg_start) {
365	3x	best_score = score;
366	3x	best_deg = d;
367	3x	best_fit.degree = d;
368	3x	best_fit.coeffs_x = basis_raw_coeffs * beta_x;
369	3x	best_fit.coeffs_y = basis_raw_coeffs * beta_y;
370	3x	best_fit.basis = basis;
371	3x	best_fit.basis_alpha = basis_alpha;
372	3x	best_fit.basis_norm2 = basis_norm2;
373	3x	best_fit.cv_error = cv_err;
374	3x	best_fit.aic = aic;
375	3x	best_fit.t = t_vec;
376	3x	best_fit.fitted = arma::join_rows(pred_x, pred_y);
377		}
378		}
379
380	3x	return best_fit;
381		}
382
383		/** @brief Compute trajectory velocities, speeds, accelerations, heading turn rates, and curvatures. */
384	1x	inline TrajectoryDerivatives eval_trajectory_derivatives(
385		const arma::vec& coeffs_x,
386		const arma::vec& coeffs_y,
387		const arma::vec& t_eval
388		) {
389	1x	size_t n = t_eval.n_elem;
390	1x	TrajectoryDerivatives d;
391	1x	d.t = t_eval;
392	1x	d.vx.set_size(n);
393	1x	d.vy.set_size(n);
394	1x	d.speed.set_size(n);
395	1x	d.ax.set_size(n);
396	1x	d.ay.set_size(n);
397	1x	d.heading_rate.set_size(n);
398	1x	d.curvature.set_size(n);
399
400	4x	for (size_t i = 0; i < n; ++i) {
401	3x	double t = t_eval[i];
402	3x	double vx = eval_poly_deriv1_scalar(coeffs_x, t);
403	3x	double vy = eval_poly_deriv1_scalar(coeffs_y, t);
404	3x	double ax = eval_poly_deriv2_scalar(coeffs_x, t);
405	3x	double ay = eval_poly_deriv2_scalar(coeffs_y, t);
406
407	3x	double spd = std::sqrt(vx * vx + vy * vy);
408	3x	double omega = vx * ay - vy * ax; // angular velocity / heading rate
409	3x	double spd3 = spd * spd * spd;
410		double kappa = (spd3 > 1e-12) ? (omega / spd3) : 0.0;
411
412	3x	d.vx[i] = vx;
413	3x	d.vy[i] = vy;
414	3x	d.speed[i] = spd;
415	3x	d.ax[i] = ax;
416	3x	d.ay[i] = ay;
417	3x	d.heading_rate[i] = omega;
418	6x	d.curvature[i] = kappa;
419		}
420
421	1x	return d;
422		}
423
424		} // namespace qe
425
426		#endif // LIBQE_TRAJECTORY_HPP