Fix structure, code style and tests in autodiff

examples/autopoint/CMakeLists.txt

@@ -11,7 +11,6 @@ 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" 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

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_functors.hpp

@@ -10,12 +10,6 @@
 
 #ifndef ODE_SOLVERS_ORACLE_FUNCTORS_HPP
 #define ODE_SOLVERS_ORACLE_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 OptimizationFunctor {
     template <
@@ -400,84 +394,4 @@ struct HessianFunctor {
 
 };
 
-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 Kernel=Cartesian<NT>;
-  using Point=typename Kernel::Point;
-  static std::function<FT(const Autopoint &, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &)> pdf;
-  // static std::function<FT(NT,NT)> f;
-  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 {
-  typedef typename Point::FT NT;
-  typedef std::vector<Point> pts;
-  typedef FunctionFunctor_internal<NT> NegativeLogprobFunctor;
-  NegativeLogprobFunctor F;
-  parameters<NT> &params;
-  GradientFunctor(parameters<NT> &params_) : params(params_){};
-  // The index i represents the state vector index
-  Point operator()(unsigned int const &i, pts 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 {
-  typedef typename Point::FT NT;
-  parameters<NT> &params;
-  FunctionFunctor(parameters<NT> &params_) : params(params_){};
-  // The index i represents the state vector index
-  typedef FunctionFunctor_internal<NT> NegativeLogprobFunctor;
-  NegativeLogprobFunctor F;
-  NT operator()(Point const &x) const {
-    NT temp = F.result(x, params.data);
-    // std::cout<<temp<<std::endl<<std::flush;
-    return temp;
-  }
-};
-};
-
 #endif