Interior point methods for SOS optimization. (#93)

* First version * fix picture * fix in doc * updated Eigen library * Bug fixes. Add option to choose floating-point type for IPM * Refactoring * Added initialization scaling * Add data typ choice for IPM * Fixes * Bug fixes and added infeasibility termination * Added README * Fixes, updated README * Fix in loading boost * Fix in doc * Edited doc * Added doc * Removed file * doc * doc * Refactoring * Add timer header file * extracted test file * Refactoring * Refactoring * added file warning * refactor * Added Copyright message * Refactoring * tmp weekly report * Many fixes. Most notably faster calculation of objective integrals via Clenshaw-Curtis algorithm. * Refactor. Remove outdated code. * Added Chebyshev Tools * Speed up initialisation via Chebyshev basis and Clenshaw Curtis algorithm * added multithreading * Refactoring. Add doc * fixes * restore old Eigen library * Adjust to new directories * delete unneccesary files * speedups * More efficient (sparse and triangular) computations. * Add sum barrier * Added Weighted SOS polynomials * Refactor barriers * Refactoring * Refactoring * Refactor barriers * Fix SDP solver * Minor * Removed spdlog. From now on user provided * Refactoring * Refactor * Refactor * Add multivariate support. WIP * Added bivariate and multivariate support * Added multithreading. Added abstraction for Product Barrier * Added MKL Support Added JSON parsing * minor * Major change. Use templated double class for all algorithms * Refactoring * Revert to original Eigen Headers. * Remove redundant Eigen Headers * Refactor * Added description * Refactoring. Fixed inheritance from Eigen/LLT * fix typo * Move .md files * fixes * Remove SDP files * Cleanup and removal of unused files * Minor. Add notes * partially addressed PR Comments * Move files * Further fixes according to PR review * * Rename Files * Add Boost Info * Remove absolute paths in CMakeLists.txt * Trigger Build
GeomScale · Feb 19, 2021 · 7cf216a · 7cf216a
1 parent f6edf61
commit 7cf216a
Show file tree

Hide file tree

Showing 37 changed files with 10,008 additions and 0 deletions.
diff --git a/examples/EnvelopeProblemSOS/CMakeLists.txt b/examples/EnvelopeProblemSOS/CMakeLists.txt
@@ -0,0 +1,86 @@
+cmake_minimum_required(VERSION 3.15)
+project(EnvelopeProblem)
+
+set(CMAKE_CXX_STANDARD 14)
+
+set(CMAKE_CXX_FLAGS_DEBUG_CUSTOM "-O0 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -DIPM_USE_DOUBLE -DIPM_DOUBLE=double")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS_DEBUG_CUSTOM}")
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -Wall -Wextra -O3 -DNDEBUG=1")
+set(CMAKE_CXX_FLAGS_RELEASE_DOUBLE_AVX "${CMAKE_CXX_FLAGS} -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -Wall -Wextra -O3 -mavx -mfma -DNDEBUG=1 -DIPM_USE_DOUBLE -DDIGITS_PRECISION=250 -DIPM_DOUBLE=double")
+
+#Environment variables MKLROOT, DYLD_LIBRARY_PATH, LIBRARY_PATH, CPATH and NLSPATH must be set to use MKL. See more information here: https://software.intel.com/content/www/us/en/develop/articles/intel-mkl-link-line-advisor.html
+set(MKL_OPTIONS "-DMKL_LP64 -m64 -I$ENV{MKLROOT}/include")
+
+set(CMAKE_CXX_FLAGS_RELEASE_DOUBLE_AVX_OPENMP "${CMAKE_CXX_FLAGS_DEBUG} -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -O3 -march=native -Xclang  -fopenmp -mavx -mfma -fno-math-errno -DNDEBUG=1 -DIPM_USE_DOUBLE -DDIGITS_PRECISION=250 -DIPM_DOUBLE=double -DEIGEN_USE_MKL_ALL=1 ${MKL_OPTIONS}")
+set(CMAKE_CXX_FLAGS_RELEASE_DOUBLE "${CMAKE_CXX_FLAGS} -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -Wall -Wextra -O3 -march=native -Xclang  -fopenmp -mavx -mfma -fno-math-errno -DNDEBUG=1 -DIPM_USE_DOUBLE -DDIGITS_PRECISION=250 -DIPM_DOUBLE=double ${MKL_OPTIONS}")
+
+file(GLOB TARGETS "*.h" "*.cpp" "../../include/sos/*.h" "../../include/sos/*.cpp" "../../include/barriers/*.h" "../../include/barriers/*.cpp"
+        "../../external/ChebTools/*"
+        "../../external/Padua/*"
+        "../../include/sos/utils.cpp"
+        )
+
+if (CMAKE_BUILD_TYPE STREQUAL "Release_double_avx_openmp")
+    set(USE_OpenMP TRUE)
+else ()
+    set(USE_OpenMP FALSE)
+endif ()
+
+message("USE_OPENMP is ${USE_OpenMP}")
+
+if (APPLE AND USE_OpenMP)
+    if (CMAKE_C_COMPILER_ID MATCHES "Clang")
+        set(OpenMP_C "${CMAKE_C_COMPILER}")
+        set(OpenMP_C_FLAGS "-fopenmp=libomp -Wno-unused-command-line-argument")
+        set(OpenMP_C_LIB_NAMES "libomp" "libgomp" "libiomp5")
+        set(OpenMP_libomp_LIBRARY ${OpenMP_C_LIB_NAMES})
+        set(OpenMP_libgomp_LIBRARY ${OpenMP_C_LIB_NAMES})
+        set(OpenMP_libiomp5_LIBRARY ${OpenMP_C_LIB_NAMES})
+    endif ()
+    if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+        set(OpenMP_CXX "${CMAKE_CXX_COMPILER}")
+        set(OpenMP_CXX_FLAGS "-Wno-unused-command-line-argument")
+        set(OpenMP_CXX_LIB_NAMES "libomp" "libgomp" "libiomp5")
+        set(OpenMP_libomp_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+        set(OpenMP_libgomp_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+        set(OpenMP_libiomp5_LIBRARY ${OpenMP_CXX_LIB_NAMES})
+        message("libomp library is in ${OpenMP_libomp_LIBRARY}")
+    endif ()
+endif ()
+
+if (USE_OpenMP)
+    find_package(OpenMP REQUIRED)
+endif (USE_OpenMP)
+
+find_package(Python2 COMPONENTS Development NumPy)
+add_executable(EnvelopeProblem ${TARGETS})
+
+if (NOT BOOST_DIR)
+    find_package(Boost 1.67.0 COMPONENTS *boost usr/local/)
+    set(BOOST_DIR ${Boost_INCLUDE_DIRS})
+    message("Boost found at ${Boost_INCLUDE_DIRS}")
+endif ()
+
+if (NOT BOOST_DIR)
+    message(FATAL_ERROR "This program requires the boost library, and will not be compiled. Set with flag -BOOST_DIR.")
+else ()
+    target_include_directories(EnvelopeProblem PRIVATE ${BOOST_DIR}
+            ../../include/sos ../../external ${SPDLOG_DIR} ../../include/sos/include ${Python2_INCLUDE_DIRS} ${Python2_NumPy_INCLUDE_DIRS}
+            ../../../Eigen/eigen
+            /usr/local/include
+            ${LLVM_INCLUDE_DIRS})
+    message("LLVM INCLUDE DIRS are ${LLVM_INCLUDE_DIRS}")
+    if (CMAKE_BUILD_TYPE STREQUAL "Release_double_avx_openmp")
+        message("Link MKL libraries")
+        target_link_options(EnvelopeProblem PRIVATE -framework Accelerate /opt/local/lib/lapack/liblapacke.dylib $ENV{MKLROOT}/lib/libmkl_intel_ilp64.a $ENV{MKLROOT}/lib/libmkl_intel_thread.a $ENV{MKLROOT}/lib/libmkl_core.a -liomp5 -lpthread -lm -ldl)
+        #Common Error on MacOs is that @rpath/libiomp5.dylib can not be linked/found. Fix is to run "install_name_tool -change @rpath/libiomp5.dylib <MKL_LIB>/lib/libiomp5.dylib ./EXECUTABLE"
+    endif ()
+    find_package(LLVM REQUIRED CONFIG)
+    message(STATUS "Found LLVM ${LLVM_PACKAGE_VERSION}")
+    message(STATUS "Using LLVMConfig.cmake in: ${LLVM_DIR}")
+    target_link_directories(EnvelopeProblem PRIVATE ${llvm_libs})
+    if (USE_OpenMP)
+        target_link_directories(EnvelopeProblem PRIVATE ${OPENMP_LIBRARIES})
+    endif ()
+    target_link_libraries(EnvelopeProblem Python2::Python Python2::NumPy)
+endif ()
diff --git a/examples/EnvelopeProblemSOS/EnvelopeProblemSOS.h b/examples/EnvelopeProblemSOS/EnvelopeProblemSOS.h
@@ -0,0 +1,107 @@
+// VolEsti (volume computation and sampling library)
+//
+// Copyright (c) 2020 Bento Natura
+//
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef EXAMPLES_ENVELOPEPROBLEMSOS_H
+#define EXAMPLES_ENVELOPEPROBLEMSOS_H
+
+//TODO: Separate the main two classes in this file.
+#include <vector>
+#include "../include/sos/NonSymmetricIPM.h"
+
+#include "matplotlib-cpp/matplotlibcpp.h"
+
+namespace plt = matplotlibcpp;
+
+template <typename IPMDouble>
+class EnvelopeProblemSOS {
+    typedef std::vector<std::pair<IPMDouble, IPMDouble> > HyperRectangle;
+    typedef Vector<IPMDouble> Vector;
+    typedef Matrix<IPMDouble> Matrix;
+    typedef Vector PolynomialSOS;
+public:
+    EnvelopeProblemSOS(unsigned num_variables, unsigned max_degree, HyperRectangle &hyperRectangle_);
+    EnvelopeProblemSOS(std::string instance_file, std::string config_file);
+
+    //FIXME: Rename as currently the degree of the polynomial remains unchanged.
+    static InterpolantVector polynomial_product(InterpolantVector p1, InterpolantVector p2) {
+        assert(p1.rows() == p2.rows());
+        InterpolantVector p = InterpolantVector::Zero(p1.rows());
+        for (Eigen::Index i = 0; i < p.rows(); ++i) {
+            for (Eigen::Index j = 0; j <= i; j++) {
+                p(i) += p1(j) * p2(i - j);
+            }
+        }
+        return p;
+    }
+
+    static InterpolantVector polynomial_product(std::vector<InterpolantVector> const poly_vec) {
+        auto len = poly_vec.size();
+        assert(not poly_vec.empty());
+        if (len == 1) {
+            return poly_vec[0];
+        }
+        if (len == 2) {
+            return polynomial_product(poly_vec[0], poly_vec[1]);
+        }
+        auto mid = len / 2;
+        auto first_it = poly_vec.begin();
+        auto mid_it = poly_vec.begin() + mid;
+        auto last_it = poly_vec.end();
+        std::vector<InterpolantVector> vec1(first_it, mid_it);
+        std::vector<InterpolantVector> vec2(mid_it, last_it);
+        return polynomial_product(polynomial_product(vec1), polynomial_product(vec2));
+    }
+
+    InterpolantVector generate_zero_polynomial();
+
+    void add_polynomial(InterpolantVector &polynomial);
+
+    Instance<IPMDouble> construct_SOS_instance();
+
+    void print_solution(Solution<IPMDouble> sol);
+
+    void plot_polynomials_and_solution(const Solution<IPMDouble> &sol);
+
+    void calculate_basis_polynomials();
+
+    InterpolantMatrix get_transformation_matrix();
+
+private:
+    unsigned _n;
+    unsigned _d;
+    unsigned _L;
+    unsigned _U;
+    InterpolantVector _objectives_vector;
+    std::vector<InterpolantVector> _polynomials_bounds;
+    HyperRectangle _hyperRectangle;
+    std::vector<InterpolantVector> _basis_polynomials;
+    std::shared_ptr<spdlog::logger> _logger;
+
+    pt::ptree _config;
+    pt::ptree _instance;
+
+    //This bool sets whether the transformation matrix from the
+    //standard monomial basis to the Lagrange basis through the
+    //Chebyshev points of second kind is calculated.
+    //Set to true if you want to save runtime but have
+    //arbitrary polynomials plotted.
+    bool _input_in_interpolant_basis = true;
+
+    bool _use_weighted_polynomials = true;
+
+    //This bool sets the degree of variables up to which the solution of the envelope problem should be plotted.
+    unsigned _plot_threshold = 100;
+
+    void compute_clenshaw_curtis_integrals();
+
+    void initialize_loggers();
+
+    void initialize_problem();
+};
+
+#include "EnvelopeProblemSOS.hpp"
+
+#endif //EXAMPLES_ENVELOPEPROBLEMSOS_H