-
Notifications
You must be signed in to change notification settings - Fork 115
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
17 changed files
with
2,832 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
Oops, something went wrong.