Enable c++17 support in tests and fix saxpy calls

examples/autopoint/CMakeLists.txt

@@ -11,8 +11,7 @@ project( VolEsti )
 
 option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON)
 option(BUILTIN_EIGEN "Use eigen from ../external" OFF)
-option(BUILTIN_AUTODIFF "Use eigen from ../external" OFF)
-option(USE_MKL "Use MKL library to build eigen" ON)
+option(USE_MKL "Use MKL library to build eigen" OFF)
 
 CMAKE_MINIMUM_REQUIRED(VERSION 3.17)
 
@@ -61,8 +60,15 @@ else()
 
 endif(DISABLE_NLP_ORACLES)
 
+option(BUILTIN_AUTODIFF "Use autodiff from ../external" OFF)
 include("../../external/cmake-files/Autodiff.cmake")
 GetAutodiff()
+if (BUILTIN_AUTODIFF)
+    include_directories (BEFORE ../../external/_deps/Autodiff)
+else ()
+    include_directories(BEFORE /usr/local/include)
+endif(BUILTIN_AUTODIFF)
+
 include("../../external/cmake-files/Eigen.cmake")
 GetEigen()
 
@@ -77,11 +83,7 @@ if (BUILTIN_EIGEN)
 else ()
     include_directories(BEFORE /usr/include/eigen3)
 endif(BUILTIN_EIGEN)
-if (BUILTIN_AUTODIFF)
-    include_directories (BEFORE ../../external/_deps/Autodiff)
-else ()
-    include_directories(BEFORE /usr/local/include)
-endif(BUILTIN_AUTODIFF)
+
 
 if (USE_MKL)
 find_library(BLAS
@@ -141,7 +143,7 @@ endif()
   add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++14 standard
   add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
   add_definitions(${CMAKE_CXX_FLAGS} "-pthread")  # optimization of the compiler
-  
+
   #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
   add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
   add_definitions(${CMAKE_CXX_FLAGS} "-DMKL_ILP64")

examples/autopoint/Gaussian_mixture_autopoint.cpp

@@ -26,6 +26,7 @@
 #include "Eigen/Eigen"
 
 #include "ode_solvers/ode_solvers.hpp"
+#include "ode_solvers/oracle_autodiff_functors.hpp"
 #include "random.hpp"
 #include "random/uniform_int.hpp"
 #include "random/normal_distribution.hpp"

examples/autopoint/MultidimensionalGaussian_mixture_autopoint.cpp

@@ -24,6 +24,7 @@
 #include "Eigen/Eigen"
 
 #include "ode_solvers/ode_solvers.hpp"
+#include "ode_solvers/oracle_autodiff_functors.hpp"
 #include "random.hpp"
 #include "random/uniform_int.hpp"
 #include "random/normal_distribution.hpp"
@@ -79,7 +80,7 @@ void run_main()
     }
     start = std::chrono::high_resolution_clock::now();
     std::cerr << (long)std::chrono::duration_cast<std::chrono::microseconds>(start - stop).count();
-    for (int i = 0; i < max_actual_draws; i++) {  
+    for (int i = 0; i < max_actual_draws; i++) {
         std::cout << hmc.x.getCoefficients().transpose() << std::endl;
         hmc.apply(rng, 3);
         samples.col(i) = hmc.x.getCoefficients();

examples/autopoint/userDefinedFunction_autopoint.cpp

@@ -22,6 +22,7 @@
 #include "Eigen/Eigen"
 
 #include "ode_solvers/ode_solvers.hpp"
+#include "ode_solvers/oracle_autodiff_functors.hpp"
 #include "random.hpp"
 #include "random/uniform_int.hpp"
 #include "random/normal_distribution.hpp"
@@ -77,7 +78,7 @@ void run_main()
         hmc.apply(rng, 3);
     }
     start = std::chrono::high_resolution_clock::now();
-    for (int i = 0; i < max_actual_draws; i++) {   
+    for (int i = 0; i < max_actual_draws; i++) {
         std::cout << hmc.x.getCoefficients().transpose() << std::endl;
         hmc.apply(rng, 3);
         samples.col(i) = hmc.x.getCoefficients();

examples/autopoint/userDefinedFunction_manual_gradient.cpp

@@ -21,6 +21,7 @@
 #include "Eigen/Eigen"
 
 #include "ode_solvers/ode_solvers.hpp"
+#include "ode_solvers/oracle_autodiff_functors.hpp"
 #include "random.hpp"
 #include "random/uniform_int.hpp"
 #include "random/normal_distribution.hpp"
@@ -130,7 +131,7 @@ void run_main() {
     hmc.apply(rng, 3);
   }
   start = std::chrono::high_resolution_clock::now();
-  for (int i = 0; i < max_actual_draws; i++) { 
+  for (int i = 0; i < max_actual_draws; i++) {
     std::cout << hmc.x.getCoefficients().transpose() << std::endl;
     hmc.apply(rng, 3);
     samples.col(i) = hmc.x.getCoefficients();

examples/correlation_matrices/CMakeLists.txt

@@ -62,6 +62,16 @@ else()
 
 endif(DISABLE_NLP_ORACLES)
 
+option(BUILTIN_AUTODIFF "Use autodiff from ../external" OFF)
+include("../../external/cmake-files/Autodiff.cmake")
+GetAutodiff()
+if (BUILTIN_AUTODIFF)
+    include_directories (BEFORE ../../external/_deps/Autodiff)
+else ()
+    include_directories(BEFORE /usr/local/include)
+endif(BUILTIN_AUTODIFF)
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") #enable the c++17 support needed by autodiff
+
 include("../../external/cmake-files/Eigen.cmake")
 GetEigen()
 
@@ -125,7 +135,6 @@ else ()
     add_compile_definitions("EIGEN_NO_DEBUG")
   endif ()
 
-  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11")  # enable C++11 standard
   add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
   #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
   add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")

examples/count-linear-extensions-using-hpolytope/CMakeLists.txt

@@ -98,8 +98,6 @@ else ()
     add_compile_definitions("EIGEN_NO_DEBUG")
   endif ()
 
-
-  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11")  # enable C++11 standard
   add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
   #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
   add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")

include/cartesian_geom/autopoint.h

@@ -1,8 +1,18 @@
-#ifndef AB9D4FD9_C418_44AC_8276_CC42540EF027
-#define AB9D4FD9_C418_44AC_8276_CC42540EF027
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2024 Vissarion Fisikopoulos
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef CARTESIAN_KERNEL_AUTOPOINT_H
+#define CARTESIAN_KERNEL_AUTOPOINT_H
+
 #include <iostream>
 #include <Eigen/Eigen>
 
+/// This class manipulates a point used for automatic differentation
+/// parameterized by a number type e.g. double
+/// \tparam T Numerical Type
 template <typename T>
 class autopoint
 {
@@ -133,7 +143,7 @@ class autopoint
         temp.coeffs=coeffs.array().log().matrix();
         return temp;
     }
-    
+
     autopoint exp() const {
         autopoint temp;
         temp.d = d;
@@ -217,20 +227,20 @@ class autopoint
     {
         return autopoint((Coeff)(coeffs * ((FT)k)));
     }
-    
+
     autopoint operator* (T k)
     {
         return autopoint((Coeff)(coeffs * ((FT)k)));
     }
     // matrix multiplication
     autopoint operator* (const autopoint& autopoint_)
     {
-        return autopoint((Coeff)(coeffs * autopoint_.getCoefficients())); 
+        return autopoint((Coeff)(coeffs * autopoint_.getCoefficients()));
     }
     // matrix multiplication with normal matrix
     autopoint operator* (const coeff& matrix_)
     {
-        return autopoint((  autopoint(matrix_) * autopoint(coeffs) )); 
+        return autopoint((  autopoint(matrix_) * autopoint(coeffs) ));
     }
 
     void operator*= (const FT k)
@@ -283,7 +293,7 @@ class autopoint
         std::cout<<"\n";
     }
 
-    static autopoint all_ones(int dim) 
+    static autopoint all_ones(int dim)
     {
       autopoint p(dim);
       for (int i = 0; i < dim; i++) p.set_coord(i, 1.0);
@@ -298,7 +308,7 @@ autopoint<K> operator* ( K k, autopoint<K> const& p)
     return p * k;
 }
 
-// matrix times autopoint 
+// matrix times autopoint
 template<typename K>
 autopoint<K> operator* (  Eigen::Matrix<K, Eigen::Dynamic,Eigen::Dynamic> matrix_, autopoint<K> const& p)
 {
@@ -308,4 +318,4 @@ autopoint<K> operator* (  Eigen::Matrix<K, Eigen::Dynamic,Eigen::Dynamic> matrix
 
 
 
-#endif /* AB9D4FD9_C418_44AC_8276_CC42540EF027 */
+#endif // CARTESIAN_KERNEL_AUTOPOINT_H

include/ode_solvers/oracle_autodiff_functors.hpp

@@ -0,0 +1,94 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP
+#define ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP
+
+#include "Eigen/Eigen"
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+#include "cartesian_geom/cartesian_kernel.h"
+#include "cartesian_geom/autopoint.h"
+
+struct AutoDiffFunctor {
+
+    template <typename NT>
+    struct parameters {
+        unsigned int order;
+        NT L; // Lipschitz constant for gradient
+        NT m; // Strong convexity constant
+        NT kappa; // Condition number
+        Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> data;
+        parameters() : order(2), L(4), m(4), kappa(1){};
+    };
+
+    template <typename NT>
+    struct FunctionFunctor_internal {
+        using Autopoint = autopoint<NT>;
+        using Coeff = typename autopoint<NT>::Coeff;
+        using FT = typename autopoint<NT>::FT;
+        using Point = typename Cartesian<NT>::Point;
+
+        static std::function<FT(const Autopoint &, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &)> pdf;
+
+        FT static result_internal(const Coeff &x, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data){
+            return pdf(x, data); //
+        }
+
+        // external interface
+        Point static differentiate(Point const &x0, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data) {
+            Autopoint x = Autopoint(x0.getCoefficients()); // cast into autopoint
+            auto x1 = x.getCoefficients();
+            Coeff y = autodiff::gradient(result_internal, autodiff::wrt(x1), autodiff::at(x1, data));
+            auto result = y.template cast<NT>();
+            return -1 * Point(result);
+        }
+
+        NT static result(Point const &x0, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data) {
+            Autopoint x = Autopoint(x0.getCoefficients()); // cast to autopoint
+            auto x1 = x.getCoefficients();
+            return result_internal(x1, data).val();
+        }
+    };
+
+    template <typename Point>
+    struct GradientFunctor {
+        using NT = typename Point::FT;
+
+        FunctionFunctor_internal<NT> F;
+        parameters<NT> &params;
+
+        GradientFunctor(parameters<NT> &params_) : params(params_){};
+
+        // The index i represents the state vector index
+        Point operator()(unsigned int const &i, std::vector<Point> const &xs, NT const &t) const {
+            // std::cout<<"calling gradient functor"<<std::flush;
+            if (i == params.order - 1) {
+                return F.differentiate(xs[0], params.data);
+            }
+            else {
+                return xs[i + 1]; // returns derivative
+            }
+        }
+    };
+
+    template <typename Point>
+    struct FunctionFunctor {
+        using NT = typename Point::FT;
+        parameters<NT> &params;
+
+        FunctionFunctor(parameters<NT> &params_) : params(params_){};
+        // The index i represents the state vector index
+        FunctionFunctor_internal<NT> F;
+
+        NT operator()(Point const &x) const {
+            return F.result(x, params.data);
+        }
+    };
+};
+
+#endif //ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP