Merge branch 'GeomScale:develop' into feature/tuning

examples/correlation_matrices/sampler.cpp

@@ -92,28 +92,6 @@ void write_to_file(std::string filename, std::vector<PointType> const& randPoint
     std::cout.rdbuf(coutbuf);
 }
 
-bool is_correlation_matrix(const MT& matrix, const double tol = 1e-8){
-    //check if all the diagonal elements are ones
-    for(int i=0 ; i<matrix.rows() ; i++)
-    {
-        if(std::abs(matrix(i, i)-1.0) > tol)
-	{
-	    return false;
-	}
-    }
-
-    //check if the matrix is positive semidefinite
-    using NT = double;
-    using MatrixType = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;
-    EigenvaluesProblems<NT, MatrixType, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;
-
-    if(solver.isPositiveSemidefinite(matrix))
-    {
-        return true;
-    }
-    return false;
-}
-
 template<typename WalkType>
 void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){
 
@@ -140,26 +118,34 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned
     std::cout << walkname << " samples uniformly "<< num_points << " correlation matrices of size " << n << " with matrix PointType" << std::endl;
     std::chrono::steady_clock::time_point start, end;
     double time;
-    std::vector<PointMT> randPoints;
+    std::list<MT> randCorMatrices;
     unsigned int walkL = 1;
 
     start = std::chrono::steady_clock::now();
 
-    uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randPoints, walkL, num_points, 0);
+    uniform_correlation_sampling_MT<WalkType, PointMT, RNGType>(n, randCorMatrices, walkL, num_points, 0);
 
     end = std::chrono::steady_clock::now();
     time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
     std::cout << "Elapsed time : " << time << " (ms)" << std::endl;
 
     int valid_points = 0;
     EigenvaluesProblems<NT, MT, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;
-    for(const auto& points : randPoints){
-        if(solver.is_correlation_matrix(points.mat)){
+    for(const auto& matrix : randCorMatrices){
+    	if(solver.is_correlation_matrix(matrix)){
             valid_points++;
-        }
+    	}
     }
+
     std::cout << "Number of valid points = " << valid_points << std::endl;
 
+    std::vector<PointMT> randPoints;
+    for(const auto &mat : randCorMatrices){
+        PointMT p;
+    	p.mat = mat;
+    	randPoints.push_back(p);
+    }
+
     write_to_file<PointMT>(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints);
 }
 

examples/crhmc_cooling_gaussians/CMakeLists.txt

@@ -0,0 +1,132 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2024 Vissarion Fisikopoulos
+# Copyright (c) 2018-2024 Apostolos Chalkis
+# Copyright (c) 2024 Vladimir Necula
+
+# Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+# Code 2024 program.
+
+# 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)
+	set(MKLROOT $ENV{HOME}/intel/mkl)
+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("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+include("../../external/cmake-files/QD.cmake")
+GetQD()
+
+# Find lpsolve library
+find_library(LP_SOLVE NAMES liblpsolve55.so 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")
+
+  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 (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 ()
+
+
+  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
+  set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
+  add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING"  ${ADDITIONAL_FLAGS})  # optimization of the compiler
+  #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
+  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(volume_example volume_example.cpp)
+  target_link_libraries(volume_example QD_LIB ${MKL_LINK} ${LP_SOLVE})
+
+endif()

examples/crhmc_cooling_gaussians/volume_example.cpp

@@ -0,0 +1,60 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
+// Copyright (c) 2024 Vladimir Necula
+
+// Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+// Code 2024 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#include "generators/known_polytope_generators.h"
+#include "random_walks/random_walks.hpp"
+#include "volume/volume_cooling_gaussians_crhmc.hpp"
+
+#include <iostream>
+#include <fstream>
+#include "misc/misc.h"
+
+typedef double NT;
+typedef Cartesian<NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt11213b, double> RandomNumberGenerator;
+typedef HPolytope<Point> HPOLYTOPE;
+
+void calculateAndVerifyVolume(HPOLYTOPE& polytope) {
+    int walk_len = 10;
+    NT e = 0.1;
+
+    RandomNumberGenerator rng(polytope.dimension());
+
+    NT volume = volume_cooling_gaussians<HPOLYTOPE, RandomNumberGenerator>(polytope, rng, e, walk_len);
+
+    std::cout << "Volume " << volume << std::endl;
+}
+
+int main() {
+
+    HPOLYTOPE simplex = generate_simplex<HPOLYTOPE>(2, false);
+    std::cout << std::endl << "Simplex: " << std::endl;
+    simplex.print();
+    calculateAndVerifyVolume(simplex);
+
+    HPOLYTOPE cube = generate_cube<HPOLYTOPE>(3, false);
+    std::cout << std::endl << "Cube: " << std::endl;
+    cube.print();
+    calculateAndVerifyVolume(cube);
+
+    HPOLYTOPE cross = generate_cross<HPOLYTOPE>(3, false);
+    std::cout << std::endl << "Cross: " << std::endl;
+    cross.print();
+    calculateAndVerifyVolume(cross);
+
+    HPOLYTOPE birkhoff = generate_birkhoff<HPOLYTOPE>(3);
+    std::cout << std::endl << "Birkhoff: " << std::endl;
+    birkhoff.print();
+    calculateAndVerifyVolume(birkhoff);
+
+    return 0;
+}

examples/sampling-hpolytope-with-billiard-walks/sampler.cpp

@@ -7,6 +7,7 @@
 #include "sampling/random_point_generators.hpp"
 #include "random_walks/random_walks.hpp"
 #include "preprocess/max_inscribed_ellipsoid.hpp"
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
 #include "convex_bodies/ellipsoid.h"
 #include "convex_bodies/hpolytope.h"
 
@@ -69,23 +70,18 @@ void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned i
     Point q(HP.dimension()); // origin
 
     // ----------- Get inscribed ellipsoid --------------------------------
-    typedef Ellipsoid<Point> EllipsoidType;
     unsigned int max_iter = 150;
     NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
     VT x0 = q.getCoefficients();
-    MT E;
     VT center;
     bool converged;
-    std::tie(E, center, converged) = max_inscribed_ellipsoid<MT>(HP.get_mat(),
-        HP.get_vec(), x0, max_iter, tol, reg);
-
-    if (!converged) // not converged
-        throw std::runtime_error("max_inscribed_ellipsoid not converged");
-
-    EllipsoidType inscribed_ellipsoid(E);
+    std::tuple<MT, VT, NT> ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::VOLUMETRIC_BARRIER>
+    (HP.get_mat(), HP.get_vec(), x0, max_iter, tol, reg);
+
+    const MT E = get<0>(ellipsoid);
     // --------------------------------------------------------------------
 
-    Generator::apply(HP, q, inscribed_ellipsoid, num_points, walk_len,
+    Generator::apply(HP, q, E, num_points, walk_len,
                      randPoints, push_back_policy, rng);
     write_to_file(filename, randPoints);
 }

include/convex_bodies/ballintersectconvex.h

@@ -50,6 +50,12 @@ class BallIntersectPolytope {
         return P.get_vec();
     }
 
+    bool is_normalized () {
+        return true;
+    }
+
+    void normalize() {}
+
     int is_in(PointType const& p) const
     {
         if (B.is_in(p)==-1)

include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp

@@ -261,6 +261,12 @@ class CorrelationSpectrahedron : public Spectrahedron<Point>{
     MT get_mat() const {
         return MT::Identity(this->d, this->d);
     }
+
+    bool is_normalized () {
+        return true;
+    }
+
+    void normalize() {}
 };
 
 #endif //VOLESTI_CONVEX_BODIES_CORRELATION_MATRICES_VOLESTI_CORRELATION_SPECTRAHEDRON_HPP