Count linear extensions via OrderPolytope volume approximation with optimized HMC oracles

R/R_Spectrahedron.R

@@ -0,0 +1,133 @@
+#' Spectrahedron Class
+#'
+#' An S4 class representing a spectrahedron defined by a Linear Matrix Inequality (LMI).
+#'
+#' @slot ptr External pointer to C++ Spectrahedron object
+#' @slot dimension Numeric: dimension of the space
+#' @slot lmi_matrices List: matrices defining the LMI (if available)
+#' @slot interior_point Numeric: coordinates of an interior point (if known)
+#'
+#' @export
+setClass(
+  "Spectrahedron",
+  slots = list(
+    ptr = "externalptr",
+    dimension = "numeric",
+    lmi_matrices = "list",
+    interior_point = "numeric"
+  ),
+  validity = function(object) {
+    if (!is.null(object@ptr) && typeof(object@ptr) != "externalptr") {
+      return("ptr must be an external pointer")
+    }
+    if (length(object@dimension) != 1 || object@dimension <= 0) {
+      return("dimension must be positive")
+    }
+    TRUE
+  }
+)
+
+#' Create Spectrahedron from SDPA Format File
+#'
+#' @param sdpa_file Character string: path to SDPA format file
+#' @param interior_point Numeric vector (optional): starting interior point
+#'
+#' @return \code{Spectrahedron} object
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#'   S <- spectrahedron_from_sdpa("path/to/file.txt")
+#' }
+spectrahedron_from_sdpa <- function(sdpa_file, interior_point = NULL) {
+
+  if (!file.exists(sdpa_file)) {
+    stop("File not found: ", sdpa_file)
+  }
+
+  ptr <- .load_spectrahedron_from_sdpa_cpp(sdpa_file)
+
+  if (is.null(ptr) || !inherits(ptr, "externalptr")) {
+    stop("Failed to load spectrahedron from file")
+  }
+
+  dim <- .get_spectrahedron_dimension_cpp(ptr)
+
+  new(
+    "Spectrahedron",
+    ptr = ptr,
+    dimension = dim,
+    lmi_matrices = list(),
+    interior_point = if (!is.null(interior_point)) as.numeric(interior_point) else numeric(0)
+  )
+}
+
+#' Create Spectrahedron from Matrices
+#'
+#' @param L0 Numeric matrix: the constant matrix L_0 in the LMI
+#' @param L_list List of numeric matrices: matrices L_1, ..., L_n
+#'
+#' @return \code{Spectrahedron} object
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#'   L0 <- matrix(c(1, 0, 0, 1), 2, 2)
+#'   L1 <- matrix(c(-1, 0, 0, 0), 2, 2)
+#'   L2 <- matrix(c(0, 0, 0, -1), 2, 2)
+#'   S <- spectrahedron_from_matrices(L0, list(L1, L2))
+#' }
+spectrahedron_from_matrices <- function(L0, L_list) {
+
+  if (!is.matrix(L0) || !is.numeric(L0)) {
+    stop("L0 must be a numeric matrix")
+  }
+
+  if (!is.list(L_list)) {
+    stop("L_list must be a list of matrices")
+  }
+
+  if (nrow(L0) != ncol(L0)) {
+    stop("L0 must be square matrix")
+  }
+
+  n_rows <- nrow(L0)
+  n_cols <- ncol(L0)
+
+  for (i in seq_along(L_list)) {
+    if (!is.matrix(L_list[[i]]) || !is.numeric(L_list[[i]])) {
+      stop("L_list[[", i, "]] must be a numeric matrix")
+    }
+    if (nrow(L_list[[i]]) != n_rows || ncol(L_list[[i]]) != n_cols) {
+      stop("All matrices must have same dimensions as L0")
+    }
+  }
+
+  n_vars <- length(L_list)
+
+  ptr <- .create_spectrahedron_from_matrices_cpp(L0, L_list)
+
+  if (is.null(ptr) || !inherits(ptr, "externalptr")) {
+    stop("Failed to create spectrahedron from matrices")
+  }
+
+  new(
+    "Spectrahedron",
+    ptr = ptr,
+    dimension = n_vars,
+    lmi_matrices = c(list(L0), L_list),
+    interior_point = numeric(0)
+  )
+}
+
+#' @export
+setMethod("dim", "Spectrahedron", function(x) {
+  as.integer(x@dimension)
+})
+
+#' @export
+setMethod("show", "Spectrahedron", function(object) {
+  cat("Spectrahedron object\n")
+  cat("  Dimension:", object@dimension, "\n")
+  cat("  Interior point:", if (length(object@interior_point) > 0) "computed" else "not set", "\n")
+})

R/R_sample_points_modification.R

@@ -0,0 +1,87 @@
+#' @rdname sample_points
+#' @description
+#'   Sample random points from a Spectrahedron using various random walk methods
+#'
+#' @param object A \code{Spectrahedron} object
+#' @param n Integer: number of points to sample
+#' @param walk_type Character: type of random walk ("RDHR", "HMC", "CDHR", "BILLIARD")
+#' @param walk_length Integer: length of each walk segment
+#' @param random_seed Integer: seed for random number generator
+#' @param starting_point Numeric vector: initial point (auto-computed if NULL)
+#' @param n_burns Integer: number of burn-in samples
+#' @param ... Additional arguments (ignored)
+#'
+#' @return Matrix: rows are sampled points (n x d matrix where d is dimension)
+#'
+#' @export
+setMethod(
+  "sample_points",
+  signature(object = "Spectrahedron"),
+  function(object,
+           n,
+           walk_type = "RDHR",
+           walk_length = 10L,
+           random_seed = 1L,
+           starting_point = NULL,
+           n_burns = 0L,
+           ...) {
+
+    if (object@dimension == 0 || is.null(object@ptr)) {
+      stop("Invalid Spectrahedron object")
+    }
+
+    n <- as.integer(n)
+    if (n <= 0L) {
+      stop("n must be a positive integer")
+    }
+
+    walk_type <- match.arg(
+      walk_type,
+      c("RDHR", "HMC", "CDHR", "BILLIARD"),
+      several.ok = FALSE
+    )
+
+    walk_length <- as.integer(walk_length)
+    if (walk_length < 1L) {
+      stop("walk_length must be >= 1")
+    }
+
+    random_seed <- as.integer(random_seed)
+    n_burns <- as.integer(n_burns)
+    if (n_burns < 0L) {
+      stop("n_burns must be >= 0")
+    }
+
+    if (!is.null(starting_point)) {
+      starting_point <- as.numeric(starting_point)
+      if (length(starting_point) != object@dimension) {
+        stop("starting_point dimension must match Spectrahedron dimension")
+      }
+    }
+
+    tryCatch({
+      sample_matrix <- sample_spectrahedra_rcpp(
+        object@ptr,
+        n = n,
+        walk_type = walk_type,
+        walk_length = walk_length,
+        seed = random_seed,
+        starting_point = starting_point,
+        n_burns = n_burns
+      )
+
+      if (!is.matrix(sample_matrix)) {
+        stop("C++ function returned invalid type")
+      }
+
+      if (nrow(sample_matrix) != object@dimension) {
+        sample_matrix <- t(sample_matrix)
+      }
+
+      return(sample_matrix)
+
+    }, error = function(e) {
+      stop("Error sampling points: ", conditionMessage(e))
+    })
+  }
+)

R/R_volume_modification.R

@@ -0,0 +1,73 @@
+#' @rdname volume
+#' @description
+#'   Compute volume of a Spectrahedron using Cooling methods
+#'
+#' @param object A \code{Spectrahedron} object
+#' @param volume_method Character: volume computation algorithm ("CB", "SOB")
+#' @param walk_length Integer: length of each random walk segment
+#' @param tolerance Numeric: relative error tolerance for volume estimate
+#' @param n_threads Integer: number of parallel threads (reserved for future)
+#' @param parameters List: additional parameters, e.g. list(seed=123)
+#' @param ... Additional arguments (ignored)
+#'
+#' @return Numeric: estimated volume of the spectrahedron
+#'
+#' @export
+setMethod(
+  "volume",
+  signature(object = "Spectrahedron"),
+  function(object, 
+           volume_method = "CB",
+           walk_length = 10L,
+           tolerance = 0.1,
+           n_threads = 1L,
+           parameters = list(),
+           ...) {
+
+    if (object@dimension == 0 || is.null(object@ptr)) {
+      stop("Invalid Spectrahedron object")
+    }
+
+    volume_method <- match.arg(
+      volume_method,
+      c("CB", "SOB"),
+      several.ok = FALSE
+    )
+
+    walk_length <- as.integer(walk_length)
+    if (walk_length < 1L) {
+      stop("walk_length must be >= 1")
+    }
+
+    if (!is.numeric(tolerance) || tolerance <= 0) {
+      stop("tolerance must be a positive number")
+    }
+
+    seed <- as.integer(parameters$seed %||% 1L)
+
+    tryCatch({
+      vol <- volume_spectrahedra_rcpp(
+        object@ptr,
+        walk_length = walk_length,
+        tolerance = tolerance,
+        seed = seed,
+        volume_method = volume_method
+      )
+
+      if (!is.numeric(vol) || !is.finite(vol)) {
+        warning("Volume computation resulted in invalid value")
+        return(NA_real_)
+      }
+
+      return(vol)
+
+    }, error = function(e) {
+      stop("Error computing volume: ", conditionMessage(e))
+    })
+  }
+)
+
+# Helper function
+`%||%` <- function(x, y) {
+  if (is.null(x)) y else x
+}

examples/count-linear-extensions/CMakeLists.txt

@@ -0,0 +1,73 @@
+# VolEsti (volume computation and sampling library)
+# Licensed under GNU LGPL.3, see LICENCE file
+
+project( VolEsti )
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+if(COMMAND cmake_policy)
+  cmake_policy(SET CMP0003 NEW)
+endif(COMMAND cmake_policy)
+
+option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON)
+option(BUILTIN_EIGEN "Use eigen from ../../external" OFF)
+
+if(DISABLE_NLP_ORACLES)
+  add_definitions(-DDISABLE_NLP_ORACLES)
+endif(DISABLE_NLP_ORACLES)
+
+include("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+# Find lpsolve library
+find_library(LP_SOLVE NAMES liblpsolve55.so liblpsolve55.a PATHS /usr/lib/lp_solve)
+
+if (NOT LP_SOLVE)
+  message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.")
+else ()
+  message(STATUS "Library lp_solve found: ${LP_SOLVE}")
+
+  set(CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+  include_directories (BEFORE ../../external)
+  include_directories (BEFORE ../../external/minimum_ellipsoid)
+  include_directories (BEFORE ../../include/generators)
+  include_directories (BEFORE ../../include/volume)
+  include_directories (BEFORE ../../include)
+  include_directories (BEFORE ../../include/convex_bodies)
+  include_directories (BEFORE ../../include/random_walks)
+  include_directories (BEFORE ../../include/annealing)
+  include_directories (BEFORE ../../include/ode_solvers)
+  include_directories (BEFORE ../../include/root_finders)
+  include_directories (BEFORE ../../include/samplers)
+  include_directories (BEFORE ../../include/lp_oracles)
+  include_directories (BEFORE ../../include/nlp_oracles)
+  include_directories (BEFORE ../../include/misc)
+  include_directories (BEFORE ../../include/optimization)
+  include_directories (BEFORE ../../include/preprocess)
+
+  # for Eigen
+  if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
+    add_compile_options(-D "EIGEN_NO_DEBUG")
+  else ()
+    add_compile_definitions("EIGEN_NO_DEBUG")
+  endif ()
+
+  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
+  add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR")
+
+  add_executable (count_le_orderpolytope count_le_orderpolytope.cpp)
+  TARGET_LINK_LIBRARIES(count_le_orderpolytope ${LP_SOLVE})
+
+endif()