Added crhmc in volesti4dingo

include/ode_solvers/implicit_midpoint.hpp

@@ -0,0 +1,180 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2022-2022 Ioannis Iakovidis
+
+// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of
+// Code 2022 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+// References
+// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with
+// Riemannian Hamiltonian
+// Monte Carlo in a Constrained Space"
+
+#ifndef IMPLICIT_MIDPOINT_HPP
+#define IMPLICIT_MIDPOINT_HPP
+#include "preprocess/crhmc/opts.h"
+#include "random_walks/crhmc/hamiltonian.hpp"
+#ifdef TIME_KEEPING
+#include <chrono>
+#endif
+
+template <typename T>
+inline std::vector<T> operator+(const std::vector<T> &v1,
+                                const std::vector<T> &v2)
+{
+  std::vector<T> result(v1.size());
+  for (int i = 0; i < v1.size(); i++) {
+    result[i] = v1[i] + v2[i];
+  }
+  return result;
+}
+template <typename T, typename Type>
+inline std::vector<T> operator*(const std::vector<T> &v, const Type alpha)
+{
+  std::vector<T> result(v.size());
+  for (int i = 0; i < v.size(); i++) {
+    result[i] = v[i] * alpha;
+  }
+  return result;
+}
+template <typename T, typename Type>
+inline std::vector<T> operator/(const std::vector<T> &v, const Type alpha)
+{
+  return v * (1 / alpha);
+}
+template <typename T>
+inline std::vector<T> operator-(const std::vector<T> &v1,
+                                const std::vector<T> &v2)
+{
+
+  return v1 + v2 * (-1.0);
+}
+template <typename Point, typename NT, typename Polytope, typename func>
+struct ImplicitMidpointODESolver {
+  using VT = typename Polytope::VT;
+  using MT = typename Polytope::MT;
+  using pts = std::vector<Point>;
+  using hamiltonian = Hamiltonian<Polytope, Point>;
+  using Opts = opts<NT>;
+
+  unsigned int dim;
+
+  NT eta;
+  int num_steps = 0;
+  NT t;
+
+  // Contains the sub-states
+  pts xs;
+  pts xs_prev;
+
+  // Function oracle
+  func F;
+  Polytope &P;
+  Opts &options;
+  VT nu;
+
+  hamiltonian ham;
+
+  bool done;
+#ifdef TIME_KEEPING
+  std::chrono::time_point<std::chrono::high_resolution_clock> start, end;
+  std::chrono::duration<double> DU_duration =
+      std::chrono::duration<double>::zero();
+  std::chrono::duration<double> approxDK_duration =
+      std::chrono::duration<double>::zero();
+#endif
+  ImplicitMidpointODESolver(NT initial_time,
+                            NT step,
+                            pts initial_state,
+                            func oracle,
+                            Polytope &boundaries,
+                            Opts &user_options) :
+                            eta(step),
+                            t(initial_time),
+                            xs(initial_state),
+                            F(oracle),
+                            P(boundaries),
+                            options(user_options),
+                            ham(hamiltonian(boundaries))
+  {
+    dim = xs[0].dimension();
+  };
+
+  void step(int k, bool accepted) {
+    pts partialDerivatives;
+#ifdef TIME_KEEPING
+    start = std::chrono::system_clock::now();
+#endif
+    partialDerivatives = ham.DU(xs);
+#ifdef TIME_KEEPING
+    end = std::chrono::system_clock::now();
+    DU_duration += end - start;
+#endif
+    xs = xs + partialDerivatives * (eta / 2);
+    xs_prev = xs;
+    done = false;
+    nu = VT::Zero(P.equations());
+    for (int i = 0; i < options.maxODEStep; i++) {
+      pts xs_old = xs;
+      pts xmid = (xs_prev + xs) / 2.0;
+#ifdef TIME_KEEPING
+      start = std::chrono::system_clock::now();
+#endif
+      partialDerivatives = ham.approxDK(xmid, nu);
+#ifdef TIME_KEEPING
+      end = std::chrono::system_clock::now();
+      approxDK_duration += end - start;
+#endif
+      xs = xs_prev + partialDerivatives * (eta);
+      NT dist = ham.x_norm(xmid, xs - xs_old) / eta;
+      NT maxdist = dist;
+      //If the estimate does not change terminate
+      if (maxdist < options.implicitTol) {
+        done = true;
+        num_steps = i;
+        break;
+      //If the estimate is very bad sample another velocity
+      } else if (maxdist > options.convergence_limit) {
+        xs = xs * std::nan("1");
+        done = true;
+        num_steps = i;
+        break;
+      }
+    }
+#ifdef TIME_KEEPING
+    start = std::chrono::system_clock::now();
+#endif
+    partialDerivatives = ham.DU(xs);
+#ifdef TIME_KEEPING
+    end = std::chrono::system_clock::now();
+    DU_duration += end - start;
+#endif
+    xs = xs + partialDerivatives * (eta / 2);
+    ham.project(xs);
+  }
+
+  void steps(int num_steps, bool accepted) {
+    for (int i = 0; i < num_steps; i++) {
+      step(i, accepted);
+    }
+  }
+
+  Point get_state(int index) { return xs[index]; }
+
+  void set_state(int index, Point p) { xs[index] = p; }
+  void print_state() {
+    for (int j = 0; j < xs.size(); j++) {
+      std::cerr << "state " << j << ": ";
+      for (unsigned int i = 0; i < xs[j].dimension(); i++) {
+        std::cerr << xs[j][i] << " ";
+      }
+      std::cerr << '\n';
+    }
+  }
+};
+
+#endif

include/ode_solvers/ode_solvers.hpp

@@ -44,6 +44,7 @@
 #include <vector>
 
 #include "ode_solvers/euler.hpp"
+#include "ode_solvers/implicit_midpoint.hpp"
 #include "ode_solvers/runge_kutta.hpp"
 #include "ode_solvers/leapfrog.hpp"
 #include "ode_solvers/richardson_extrapolation.hpp"

include/preprocess/crhmc/analytic_center.h

@@ -0,0 +1,178 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2022-2022 Ioannis Iakovidis
+
+// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of
+// Code 2022 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+// References
+// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with
+// Riemannian Hamiltonian
+// Monte Carlo in a Constrained Space"
+#ifndef ANALYTIC_CENTER_H
+#define ANALYTIC_CENTER_H
+#include "Eigen/Eigen"
+#include "PackedCSparse/PackedChol.h"
+#include "preprocess/crhmc/crhmc_utils.h"
+#include "preprocess/crhmc/opts.h"
+#include <fstream>
+#include <iostream>
+#include <vector>
+#ifndef SIMD_LEN
+#define SIMD_LEN 0
+#endif
+const size_t chol_k2 = (SIMD_LEN == 0) ? 1 : SIMD_LEN;
+
+using NT = double;
+using MT = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;
+using VT = Eigen::Matrix<NT, Eigen::Dynamic, 1>;
+using SpMat = Eigen::SparseMatrix<NT>;
+using CholObj = PackedChol<chol_k2, int>;
+using Triple = Eigen::Triplet<double>;
+using Tx = FloatArray<double, chol_k2>;
+using Opts = opts<NT>;
+/*This function computes the analytic center of the polytope*/
+//And detects additional constraint that need to be added
+// x - It outputs the minimizer of min f(x) subjects to {Ax=b}
+// C - detected constraint matrix
+//     If the domain ({Ax=b} intersect dom(f)) is not full dimensional in {Ax=b}
+//     because of the dom(f), the algorithm will detect the collapsed dimension
+//     and output the detected constraint C x = d
+// d - detected constraint vector
+template <typename Polytope>
+std::tuple<VT, SpMat, VT> analytic_center(SpMat const &A, VT const &b, Polytope &f, Opts const &options, VT x = VT::Zero(0, 1))
+{
+  // initial conditions
+  int n = A.cols();
+  int m = A.rows();
+  if (x.rows() == 0 || !f.barrier.feasible(x))
+  {
+    x = f.barrier.center;
+  }
+
+  VT lambda = VT::Zero(n, 1);
+  int fullStep = 0;
+  NT tConst = 0;
+  NT primalErr = std::numeric_limits<NT>::max();
+  NT dualErr = std::numeric_limits<NT>::max();
+  NT primalErrMin = std::numeric_limits<NT>::max();
+  NT primalFactor = 1;
+  NT dualFactor = 1 + b.norm();
+  std::vector<int> idx;
+
+  CholObj solver = CholObj(transform_format<SpMat,NT,int>(A));
+
+  for (int iter = 0; iter < options.ipmMaxIter; iter++)
+  {
+    std::pair<VT, VT> pair_analytic_oracle = f.analytic_center_oracle(x);
+    VT grad = pair_analytic_oracle.first;
+    VT hess = pair_analytic_oracle.second;
+
+    // compute the residual
+    VT rx = lambda - grad;
+    VT rs = b - A * x;
+
+    // check stagnation
+    primalErrMin = std::min(primalErr, primalErrMin);
+    primalErr = rx.norm() / primalFactor;
+    NT dualErrLast = dualErr;
+    dualErr = rs.norm() / dualFactor;
+    bool feasible = f.barrier.feasible(x);
+    //Compare the dual and primal error to the last and minimum so far
+    if ((dualErr > (1 - 0.9 * tConst) * dualErrLast) ||
+        (primalErr > 10 * primalErrMin) || !feasible)
+    {
+      VT dist = f.barrier.boundary_distance(x);
+      NT th = options.ipmDistanceTol;
+      visit_lambda(dist, [&idx, th](double v, int i, int j)
+                   {
+        if (v < th)
+          idx.push_back(i); });
+      if (idx.size() > 0)
+      {
+        break;
+      }
+    }
+
+    // compute the step direction
+    VT Hinv = hess.cwiseInverse();
+    solver.decompose((Tx *)Hinv.data());
+    VT out(m, 1);
+    solver.solve((Tx *)rs.data(), (Tx *)out.data());
+    VT dr1 = A.transpose() * out;
+    VT in = A * Hinv.cwiseProduct(rx);
+    solver.solve((Tx *)in.data(), (Tx *)out.data());
+
+    VT dr2 = A.transpose() * out;
+    VT dx1 = Hinv.cwiseProduct(dr1);
+    VT dx2 = Hinv.cwiseProduct(rx - dr2);
+
+    // compute the step size
+    VT dx = dx1 + dx2;
+    NT tGrad = std::min(f.barrier.step_size(x, dx), 1.0);
+    dx = dx1 + tGrad * dx2;
+    NT tConst = std::min(0.99 * f.barrier.step_size(x, dx), 1.0);
+    tGrad = tGrad * tConst;
+
+    // make the step
+    x = x + tConst * dx;
+    lambda = lambda - dr2;
+
+    if (!f.barrier.feasible(x))
+    {
+      break;
+    }
+    //If we have have converged
+    if (tGrad == 1)
+    {
+      //do some more fullStep
+      fullStep = fullStep + 1;
+      if (fullStep > log(dualErr / options.ipmDualTol) && fullStep > options.min_convergence_steps)
+      {
+        break;
+      }
+    }
+    else
+    {
+      fullStep = 0;
+    }
+  }
+  SpMat C;
+  VT d;
+  if (idx.size() == 0)
+  {
+    VT dist = f.barrier.boundary_distance(x);
+    NT th = options.ipmDistanceTol;
+    visit_lambda(dist, [&idx, th](double v, int i, int j)
+                 {
+      if (v < th)
+        idx.push_back(i); });
+  }
+
+  if (idx.size() > 0)
+  {
+    C.resize(idx.size(), n);
+    std::pair<VT, VT> pboundary = f.barrier.boundary(x);
+    VT A_ = pboundary.first;
+    VT b_ = pboundary.second;
+    A_ = A_(idx);
+    std::vector<Triple> sparseIdx;
+    for (int i = 0; i < idx.size(); i++)
+    {
+      sparseIdx.push_back(Triple(i, i, A_(i)));
+    }
+    C.setFromTriplets(sparseIdx.begin(), sparseIdx.end());
+    d = b_(idx);
+  }
+  else
+  {
+    C = MT::Zero(0, n).sparseView();
+    d = VT::Zero(0, 1);
+  }
+  return std::make_tuple(x, C, d);
+}
+#endif