ctmc_memory_gate.h

Go to the documentation of this file.

#ifndef LINE_API_MC_CTMC_MEMORY_GATE_H
#define LINE_API_MC_CTMC_MEMORY_GATE_H
 
/**
 * @file
 * @ingroup api_mc
 * Host-aware memory pre-gate for SolverCTMC.
 *
 * Port of `matlab/src/api/mc/ctmc_memory_gate.m`, mirrored by the JAR
 * (`MemoryGuard.gate`) and native Python (`ctmc_memory_gate`).
 *
 * THE POLICY IS REFUSAL, NOT A WARNING. When the predicted peak exceeds the
 * safe budget the solve is REFUSED with an error naming the two numbers, and
 * only an explicit `force` downgrades that to a warning. A warning the caller
 * can walk past is the wrong contract here: the failure mode it precedes is the
 * OOM killer taking the whole process, which produces no diagnostic at all and,
 * observed on 2026-08-01, took down five sessions and the terminal hosting
 * them. Refusing costs the caller a re-run with another solver; not refusing
 * costs it everything running on the machine.
 *
 * MDD IS NEVER GATED. It is served by its own analyzer and returns before any
 * explicit-generator path, so the gate is not on its route.
 * That exemption is load-bearing rather than incidental: `mdd` holds the
 * reachable set in a decision diagram whose size is governed by the diagram's
 * compression, not by the state count, so it routinely solves models this
 * estimator scores at exp(200). Gating on explicit size would clamp exactly the
 * method that exists to beat it. Keep any new non-explicit method on the same
 * side of the gate.
 *
 * The predictor is the power law of the reference, `bytes = alpha * N^beta`,
 * evaluated in log space. As in the reference, alpha and beta are CALIBRATED
 * once per host (`ctmc_get_calibration`): the sparse LU of two lattice
 * generators (40x40 and 80x80) is profiled for its fill-in bytes and its time,
 * a power law is fitted through the two points, and the fit is cached in a
 * JSON file in the temporary directory keyed by a machine signature. When the
 * profile cannot run the reference's FALLBACK coefficients (alpha = 128,
 * beta = 1.3) are used instead, which is the same branch the reference takes
 * then. The factorization here is a banded sparse LU in natural order, which
 * is deterministic across hosts; only the timing coefficients depend on the
 * machine, and the gate does not read those.
 */
 
#include <chrono>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <mutex>
#include <sstream>
#include <string>
#include <thread>
#include <vector>
 
#if defined(__unix__) || defined(__APPLE__)
#include <sys/utsname.h>
#include <unistd.h>
#endif
 
#include "json.hpp"
#include "line/util/error.h"
 
namespace line {
namespace mc {
 
/// 8 bytes of value plus 8 amortized for the index, per stored nonzero.
constexpr double CTMC_BYTES_PER_NZ = 16.0;
/// Fraction of available memory the solver may target.
constexpr double CTMC_DEFAULT_SAFETY_FRACTION = 0.6;
/// Fallback power-law coefficients, identical to MATLAB and Python.
constexpr double CTMC_FALLBACK_ALPHA = CTMC_BYTES_PER_NZ * 8.0;
constexpr double CTMC_FALLBACK_BETA = 1.3;
/// Conservative available-memory default when the host probe fails.
constexpr double CTMC_FALLBACK_AVAIL_BYTES = 1.0 * 1024.0 * 1024.0 * 1024.0;
/// The two lattice sides the reference profiles (G1, G2 of ctmc_get_calibration).
constexpr int CTMC_CALIB_G1 = 40;
constexpr int CTMC_CALIB_G2 = 80;
 
/** The gate verdict, plus the message the caller reports either way. */
struct CtmcGateResult {
    bool ok = true;
    std::string message;
    double predicted_gb = 0.0;
    double budget_gb = 0.0;
};
 
/** The per-host power-law fit, the reference's `calib` struct. */
struct CtmcCalibration {
    std::string sig;
    double alpha_mem = CTMC_FALLBACK_ALPHA;
    double beta_mem = CTMC_FALLBACK_BETA;
    double alpha_t = 0.0;
    double beta_t = 1.0;
    double timestamp = 0.0;
};
 
/**
 * Available physical memory in bytes.
 *
 * Never throws: a probe failure returns the conservative fallback, because a
 * gate that cannot measure the host must not thereby become permissive.
 */
inline double ctmc_available_memory_bytes() {
#if defined(_SC_AVPHYS_PAGES) && defined(_SC_PAGE_SIZE)
    const long pages = ::sysconf(_SC_AVPHYS_PAGES);
    const long page_size = ::sysconf(_SC_PAGE_SIZE);
    if (pages > 0 && page_size > 0)
        return static_cast<double>(pages) * static_cast<double>(page_size);
#endif
    return CTMC_FALLBACK_AVAIL_BYTES;
}
 
namespace gate_detail {
 
/** The reference's `ctmc_machine_signature`: arch | cores | platform. */
inline std::string machine_signature() {
    std::string arch =
#if defined(__x86_64__) || defined(_M_X64)
        "x86_64";
#elif defined(__aarch64__) || defined(_M_ARM64)
        "arm64";
#elif defined(__i386__) || defined(_M_IX86)
        "x86";
#else
        "unknown";
#endif
    unsigned ncores = std::thread::hardware_concurrency();
    if (ncores == 0) ncores = 1;
    std::string platform = "unknown";
#if defined(__unix__) || defined(__APPLE__)
    struct utsname u;
    if (::uname(&u) == 0) platform = std::string(u.sysname) + " " + u.machine;
#elif defined(_WIN32)
    platform = "Windows";
#endif
    return arch + "|" + std::to_string(ncores) + "|" + platform;
}
 
/** The reference's cache file, `fullfile(tempdir, 'line_ctmc_calib_<lang>.json')`. */
inline std::string cache_file() {
    const char* base = std::getenv("TMPDIR");
    std::string dir = (base && *base) ? std::string(base) : std::string("/tmp");
    if (!dir.empty() && dir.back() == '/') dir.pop_back();
    return dir + "/line_ctmc_calib_cpp.json";
}
 
/**
 * The reference's `ctmc_profile_point`: sparse LU of the (g^2 - 1)-state
 * lattice generator, returning n, the fill-in bytes and the elapsed seconds.
 *
 * The lattice is the 4-neighbour grid in column-major numbering, so its
 * bandwidth is g and the no-pivot banded LU keeps every fill-in within the
 * band; nnz(L) counts the unit diagonal, as MATLAB's `lu` does.
 */
inline void profile_point(int g, double bytes_per_nz, std::size_t& n, double& bytes,
                          double& secs) {
    const std::size_t ntot = static_cast<std::size_t>(g) * static_cast<std::size_t>(g);
    const std::size_t N = ntot - 1;  // drop the last state -> nonsingular
    n = N;
    const std::size_t bw = static_cast<std::size_t>(g);
    const std::size_t w = bw + 1;
    std::vector<double> U(N * w, 0.0);  // U[i*w + (j - i)], j in [i, i + bw]
    std::vector<double> work(2 * bw + 1, 0.0);
    std::size_t nnzL = 0, nnzU = 0;
    const std::chrono::steady_clock::time_point t0 = std::chrono::steady_clock::now();
    for (std::size_t i = 0; i < N; ++i) {
        const std::size_t lo = i >= bw ? i - bw : 0;
        const std::size_t hi = std::min(N - 1, i + bw);
        std::fill(work.begin(), work.end(), 0.0);
        // row i of A = Q(1:n-1, 1:n-1): the full-lattice degree on the diagonal,
        // 1 on each neighbour still inside the kept block
        const std::size_t r = i % bw, c = i / bw;
        int deg = 0;
        if (c > 0) ++deg;
        if (c + 1 < bw) ++deg;
        if (r > 0) ++deg;
        if (r + 1 < bw) ++deg;
        work[i - lo] = -static_cast<double>(deg);
        if (c > 0) work[(i - bw) - lo] = 1.0;
        if (c + 1 < bw && i + bw <= hi) work[(i + bw) - lo] = 1.0;
        if (r > 0) work[(i - 1) - lo] = 1.0;
        if (r + 1 < bw && i + 1 <= hi) work[(i + 1) - lo] = 1.0;
        for (std::size_t k = lo; k < i; ++k) {
            const double a = work[k - lo];
            if (a == 0.0) continue;
            const double ukk = U[k * w];
            const double lik = a / ukk;
            ++nnzL;
            const std::size_t kh = std::min(N - 1, k + bw);
            for (std::size_t j = k + 1; j <= kh; ++j) work[j - lo] -= lik * U[k * w + (j - k)];
        }
        ++nnzL;  // the unit diagonal of L
        for (std::size_t j = i; j <= hi; ++j) {
            const double v = work[j - lo];
            U[i * w + (j - i)] = v;
            if (v != 0.0) ++nnzU;
        }
    }
    const std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();
    secs = std::chrono::duration<double>(t1 - t0).count();
    bytes = bytes_per_nz * static_cast<double>(nnzL + nnzU);
}
 
/** The reference's `ctmc_fit_power_law`; false on a degenerate pair. */
inline bool fit_power_law(double x1, double y1, double x2, double y2, double& alpha,
                          double& beta) {
    if (x1 <= 0 || x2 <= 0 || y1 <= 0 || y2 <= 0 || x1 == x2) return false;
    beta = std::log(y2 / y1) / std::log(x2 / x1);
    alpha = y1 / std::pow(x1, beta);
    return true;
}
 
inline double now_posix() {
    return std::chrono::duration<double>(std::chrono::system_clock::now().time_since_epoch())
        .count();
}
 
}  // namespace gate_detail
 
/**
 * The reference's `ctmc_get_calibration`: the cached per-host fit, profiled
 * once when the cache is missing or belongs to another machine, and the
 * fallback coefficients when the profile fails.
 */
inline const CtmcCalibration& ctmc_get_calibration(bool verbose = false) {
    static CtmcCalibration calib;
    static bool ready = false;
    static std::mutex mu;
    std::lock_guard<std::mutex> lock(mu);
    if (ready) return calib;
    const std::string sig = gate_detail::machine_signature();
    const std::string cachefile = gate_detail::cache_file();
 
    // the cached fit, when it is this machine's
    try {
        std::ifstream in(cachefile);
        if (in) {
            nlohmann::json data;
            in >> data;
            if (data.contains("sig") && data["sig"].is_string() && data["sig"] == sig &&
                data.contains("alpha_mem") && data.contains("beta_mem")) {
                calib.sig = sig;
                calib.alpha_mem = data["alpha_mem"].get<double>();
                calib.beta_mem = data["beta_mem"].get<double>();
                calib.alpha_t = data.value("alpha_t", 0.0);
                calib.beta_t = data.value("beta_t", 1.0);
                calib.timestamp = data.value("timestamp", 0.0);
                ready = true;
                return calib;
            }
        }
    } catch (...) {
    }
 
    bool fitted = false;
    try {
        std::size_t n1 = 0, n2 = 0;
        double b1 = 0, b2 = 0, t1 = 0, t2 = 0, am = 0, bm = 0, at = 0, bt = 0;
        gate_detail::profile_point(CTMC_CALIB_G1, CTMC_BYTES_PER_NZ, n1, b1, t1);
        gate_detail::profile_point(CTMC_CALIB_G2, CTMC_BYTES_PER_NZ, n2, b2, t2);
        const bool okm = gate_detail::fit_power_law(static_cast<double>(n1), b1,
                                                    static_cast<double>(n2), b2, am, bm);
        const bool okt = gate_detail::fit_power_law(static_cast<double>(n1), std::max(t1, 1e-9),
                                                    static_cast<double>(n2), std::max(t2, 1e-9),
                                                    at, bt);
        if (!okm || !okt) throw NumericError(ErrorKind{"Degenerate"}, "degenerate power-law fit");
        calib.sig = sig;
        calib.alpha_mem = am;
        calib.beta_mem = bm;
        calib.alpha_t = at;
        calib.beta_t = bt;
        calib.timestamp = gate_detail::now_posix();
        fitted = true;
        try {
            nlohmann::json data;
            data["sig"] = calib.sig;
            data["alpha_mem"] = calib.alpha_mem;
            data["beta_mem"] = calib.beta_mem;
            data["alpha_t"] = calib.alpha_t;
            data["beta_t"] = calib.beta_t;
            data["timestamp"] = calib.timestamp;
            std::ofstream out(cachefile);
            if (out) out << data.dump();
        } catch (...) {
        }
        if (verbose) std::printf("CTMC calibration: bytes ~ %.3g*N^%.3f\n", am, bm);
    } catch (const std::exception& e) {
        if (verbose)
            std::printf("CTMC calibration failed (%s); using fallback model\n", e.what());
    } catch (...) {
        if (verbose) std::printf("CTMC calibration failed; using fallback model\n");
    }
    if (!fitted) {
        calib.sig = sig;
        calib.alpha_mem = CTMC_FALLBACK_ALPHA;
        calib.beta_mem = CTMC_FALLBACK_BETA;
        calib.alpha_t = 0.0;
        calib.beta_t = 1.0;
        calib.timestamp = gate_detail::now_posix();
    }
    ready = true;
    return calib;
}
 
/**
 * Decide whether a state space of log-size `log_nstates` can be solved here.
 *
 * @param log_nstates    natural log of the worst-case state count
 * @param force          true downgrades a refusal to a warning
 * @param safety_fraction fraction of available memory the solve may target
 * @param verbose        print the calibration and the forced / near-budget notes,
 *                       the reference's third argument (last here so the
 *                       existing three-argument callers keep their meaning)
 */
inline CtmcGateResult ctmc_memory_gate(double log_nstates, bool force = false,
                                       double safety_fraction = CTMC_DEFAULT_SAFETY_FRACTION,
                                       bool verbose = false) {
    CtmcGateResult res;
    const double avail = ctmc_available_memory_bytes();
    const double budget = safety_fraction * avail;
    const CtmcCalibration& calib = ctmc_get_calibration(verbose);
 
    // Compared in LOG space: the predicted byte count of an intractable model
    // overflows a double long before it exceeds the budget, and inf > x is a
    // comparison that has already lost the margin it was meant to report.
    const double log_pred = std::log(calib.alpha_mem) + calib.beta_mem * log_nstates;
    const double log_budget = std::log(budget > 1.0 ? budget : 1.0);
    const double gb = 1024.0 * 1024.0 * 1024.0;
 
    res.predicted_gb = std::exp(log_pred < 700.0 ? log_pred : 700.0) / gb;
    res.budget_gb = budget / gb;
    if (log_pred > log_budget) {
        std::ostringstream os;
        os.setf(std::ios::fixed);
        os.precision(2);
        os << "CTMC predicted peak memory ~" << res.predicted_gb
           << " GB exceeds the safe budget ~" << res.budget_gb << " GB ("
           << static_cast<long long>(std::floor(100.0 * safety_fraction + 0.5)) << "% of "
           << avail / gb
           << " GB available). Reduce the state space (e.g. lower 'cutoff'), use the 'mdd' method, "
              "use another solver (MVA/NC/FLD), or set force=true to override.";
        res.message = os.str();
        if (!force) {
            res.ok = false;
            return res;
        }
        if (verbose) std::printf("Warning (forced): %s\n", res.message.c_str());
    } else if (verbose && log_pred > std::log(std::max(0.5 * budget, 1.0))) {
        std::printf("CTMC predicted peak memory ~%.2f GB (budget ~%.2f GB).\n", res.predicted_gb,
                    res.budget_gb);
    }
    return res;
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_MEMORY_GATE_H