Full dimensional polytope computation

examples/full_dimensional_polytope/CMakeLists.txt

@@ -0,0 +1,111 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2025 Vissarion Fisikopoulos
+# Copyright (c) 2018-2025 Apostolos Chalkis
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+project( VolEsti )
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+# Locate Intel MKL root (in case it is enabled)
+if (APPLE)
+	set(MKLROOT /opt/intel/oneapi/mkl/latest)
+elseif(UNIX)
+  set(MKLROOT /opt/intel/oneapi/mkl/latest)
+endif()
+
+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)
+option(USE_MKL "Use MKL library to build eigen" OFF)
+
+if(DISABLE_NLP_ORACLES)
+  add_definitions(-DDISABLE_NLP_ORACLES)
+else()
+  find_library(IFOPT NAMES libifopt_core.so PATHS /usr/local/lib)
+  find_library(IFOPT_IPOPT NAMES libifopt_ipopt.so PATHS /usr/local/lib)
+  find_library(GMP NAMES libgmp.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(MPSOLVE NAMES libmps.so PATHS /usr/local/lib)
+  find_library(PTHREAD NAMES libpthread.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(FFTW3 NAMES libfftw3.so.3 PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+
+  if (NOT IFOPT)
+    message(FATAL_ERROR "This program requires the ifopt library, and will not be compiled.")
+  elseif (NOT GMP)
+    message(FATAL_ERROR "This program requires the gmp library, and will not be compiled.")
+  elseif (NOT MPSOLVE)
+    message(FATAL_ERROR "This program requires the mpsolve library, and will not be compiled.")
+  elseif (NOT FFTW3)
+    message(FATAL_ERROR "This program requires the fftw3 library, and will not be compiled.")
+  else()
+    message(STATUS "Library ifopt found: ${IFOPT}")
+    message(STATUS "Library gmp found: ${GMP}")
+    message(STATUS "Library mpsolve found: ${MPSOLVE}")
+    message(STATUS "Library fftw3 found:" ${FFTW3})
+  endif(NOT IFOPT)
+endif(DISABLE_NLP_ORACLES)
+
+# Include CMake modules for dependencies
+include("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+include("../../external/cmake-files/SuiteSparse.cmake")
+GetSuiteSparse()
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+if (USE_MKL)
+  find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib)
+  find_library(GFORTRAN NAME libgfortran.dylib PATHS /usr/local/Cellar/gcc/10.2.0_4/lib/gcc/10)
+  find_library(LAPACK NAME liblapack.dylib PATHS /usr/lib)
+  find_library(OPENMP NAME libiomp5.dylib PATHS /opt/intel/oneapi/compiler/2021.1.1/mac/compiler/lib)
+
+  include_directories (BEFORE ${MKLROOT}/include)
+  set(PROJECT_LIBS ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${GFORTRAN_LIBRARIES})
+  set(MKL_LINK "-L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl")
+  add_definitions(-DEIGEN_USE_MKL_ALL)
+else()
+  set(MKL_LINK "")
+endif(USE_MKL)
+
+# Include directories
+include_directories (BEFORE ../../external)
+include_directories (BEFORE ../../external/minimum_ellipsoid)
+include_directories (BEFORE ../../include/)
+
+# 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 ()
+
+# Compiler flags
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
+add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization
+add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
+add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")
+add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR")
+
+# Build the example executable
+add_executable(full_dimensional_example full_dimensional_example.cpp)
+
+# Link libraries (including SuiteSparse for the full_dimensional_polytope function)
+target_link_libraries(full_dimensional_example 
+    lp_solve 
+    ${SUITESPARSE_LIBRARIES} 
+    ${MKL_LINK}
+)

examples/full_dimensional_polytope/README.md

@@ -0,0 +1,79 @@
+# Full Dimensional Polytope Examples
+
+This directory contains examples demonstrating the `compute_full_dimensional_polytope` function from `include/preprocess/full_dimensional_polytope.hpp`.
+
+## Overview
+
+The `compute_full_dimensional_polytope` function transforms a potentially lower-dimensional polytope defined by:
+- **Equality constraints**: `Aeq * x = beq`
+- **Inequality constraints**: `A * x <= b`
+
+Into a full-dimensional polytope: `A_full * y <= b_full`
+
+The transformation provides:
+- **shift**: A point satisfying the equality constraints
+- **N**: Nullspace basis matrix with orthonormal columns
+- **Mapping**: `x = shift + N * y` (from reduced space to original space)
+
+## Building the Example
+
+From this directory:
+
+```bash
+mkdir build
+cd build
+cmake ..
+make
+```
+
+## Running the Example
+
+```bash
+./full_dimensional_example
+```
+
+The program will run through 6 examples, pausing between each for you to review the output.
+
+## Output Format
+
+For each example, the program displays:
+
+1. **Description**: What the example demonstrates
+2. **Input**: The equality and inequality constraints
+3. **Results**:
+   - Original and reduced dimensions
+   - Shift vector
+   - Nullspace basis N
+   - Output polytope dimensions
+   - Orthonormality verification
+
+## Key Features Demonstrated
+
+- ✅ Automatic dimension reduction
+- ✅ Orthonormal nullspace basis
+- ✅ Handling of sparse constraints
+- ✅ Multiple equality constraints
+- ✅ Extreme reductions (nD → 1D)
+- ✅ Real-world polytopes (Birkhoff)
+
+## Mathematical Background
+
+Given a polytope `P = {x : Aeq*x = beq, A*x <= b}` which is lower-dimensional (rank(Aeq) > 0), the function computes:
+
+1. **Shift vector**: Solves `Aeq * shift = beq` using SPQR
+2. **Nullspace basis**: Computes `N` such that `Aeq * N = 0` via QR factorization
+3. **Transformed polytope**: `P' = {y : A*N*y <= b - A*shift}`
+
+Any point `y` in `P'` maps to a feasible point in `P` via: `x = shift + N*y`
+
+## Dependencies
+
+- Eigen3 (matrix operations)
+- SuiteSparse (SPQR for QR factorization, CHOLMOD for linear systems)
+- Boost (for polytope generators)
+- lp_solve (for polytope operations)
+
+## See Also
+
+- **Header file**: `include/preprocess/full_dimensional_polytope.hpp`
+- **Tests**: `test/full_dimensional_polytope_test.cpp`