Revision 96344048632acd7871b0a92ac903fc07575517f4 authored by TolisChal on 07 December 2020, 16:23:16 UTC, committed by TolisChal on 07 December 2020, 16:23:16 UTC
VpolyCG_test.cpp
// VolEsti (volume computation and sampling library)
// Copyright (c) 2012-2018 Vissarion Fisikopoulos
// Copyright (c) 2018 Apostolos Chalkis
// Licensed under GNU LGPL.3, see LICENCE file
#include "doctest.h"
#include <unistd.h>
#include "Eigen/Eigen"
#include "random.hpp"
#include "random/uniform_int.hpp"
#include "random/normal_distribution.hpp"
#include "random/uniform_real_distribution.hpp"
#include "volume.h"
#include "known_polytope_generators.h"
#include <string>
#include <typeinfo>
template <typename NT>
NT factorial(NT n)
{
return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
}
template <typename NT, class RNGType, class Polytope>
void test_CV_volume(Polytope &VP, NT expected, NT tolerance=0.25)
{
typedef typename Polytope::PolytopePoint Point;
// Setup the parameters
int n = VP.dimension();
int walk_len=1;
int nexp=1, n_threads=1;
NT e=0.2, err=0.0000000001;
NT C=2.0,ratio,frac=0.1,delta=-1.0;
int rnum = std::pow(e,-2) * 400 * n * std::log(n);
int N = 500 * ((int) C) + ((int) (n * n / 2));
int W = 4*n*n+500;
ratio = 1.0-1.0/(NT(n));
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
RNGType rng(seed);
boost::normal_distribution<> rdist(0,1);
boost::random::uniform_real_distribution<>(urdist);
boost::random::uniform_real_distribution<> urdist1(-1,1);
//compute the chebychev ball
std::pair<Point,NT> CheBall;
// Estimate the volume
std::cout << "Number type: " << typeid(NT).name() << std::endl;
NT vol = 0;
unsigned int const num_of_exp = 10;
CheBall = VP.ComputeInnerBall();
Polytope VP2 = VP;
for (unsigned int i=0; i<num_of_exp; i++)
{
VP.init(n,VP2.get_mat(),VP2.get_vec());
vars<NT, RNGType> var2(rnum,n,10 + n/10,n_threads,err,e,0,0,0,0,0.0,rng,
urdist,urdist1,-1.0,false,false,false,false,false,false,true,false,false);
vars_g<NT, RNGType> var1(n,walk_len,N,W,1,e,CheBall.second,rng,C,frac,ratio,delta,
false,false,false,false,false,false,true,false);
vol += volume_gaussian_annealing(VP, var1, var2, CheBall);
}
NT error = std::abs(((vol/num_of_exp)-expected))/expected;
std::cout << "Computed volume (average) = " << vol/num_of_exp << std::endl;
std::cout << "Expected volume = " << expected << std::endl;
CHECK(error < tolerance);
}
template <typename NT>
void call_test_cube(){
typedef Cartesian<NT> Kernel;
typedef typename Kernel::Point Point;
typedef boost::mt19937 RNGType;
typedef VPolytope<Point, RNGType > Vpolytope;
Vpolytope P;
std::cout << "--- Testing volume of V-cube3" << std::endl;
P = generate_cube<Vpolytope>(3, true);
test_CV_volume<NT, RNGType>(P, 8.0);
std::cout << "--- Testing volume of V-cube4" << std::endl;
P = generate_cube<Vpolytope>(4, true);
test_CV_volume<NT, RNGType>(P, 16.0);
}
template <typename NT>
void call_test_cross(){
typedef Cartesian<NT> Kernel;
typedef typename Kernel::Point Point;
typedef boost::mt19937 RNGType;
typedef VPolytope<Point, RNGType > Vpolytope;
Vpolytope P;
std::cout << "--- Testing volume of V-cross4" << std::endl;
P = gen_cross<Vpolytope>(4, true);
test_CV_volume<NT, RNGType>(P, 0.6666667);
//if(typeid(NT)== typeid(double)) {
// std::cout << "--- Testing volume of V-cross5" << std::endl;
// P = gen_cross<Vpolytope>(5, true);
// test_CV_volume<NT, RNGType>(P, 0.26666667);
// }
}
template <typename NT>
void call_test_simplex() {
typedef Cartesian<NT> Kernel;
typedef typename Kernel::Point Point;
typedef boost::mt19937 RNGType;
typedef VPolytope<Point, RNGType > Vpolytope;
Vpolytope P;
std::cout << "--- Testing volume of V-simplex5" << std::endl;
P = gen_simplex<Vpolytope>(5, true);
test_CV_volume<NT, RNGType>(P, 1.0 / factorial(5.0));
//if(typeid(NT)== typeid(double)) {
//std::cout << "--- Testing volume of V-simplex7" << std::endl;
//P = gen_simplex<Vpolytope>(7, true);
//test_CV_volume<NT, RNGType>(P, 1.0 / factorial(7.0));
//}
}
TEST_CASE("cube") {
call_test_cube<double>();
//call_test_cube<float>();
//call_test_cube<long double>();
}
TEST_CASE("cross") {
call_test_cross<double>();
//call_test_cross<float>();
// call_test_cross<long double>();
}
TEST_CASE("simplex") {
call_test_simplex<double>();
//call_test_simplex<float>();
//call_test_simplex<long double>();
}
Computing file changes ...