Skip to content

Commit

Permalink
New rounding implementation ans volume approximation for V-poytopes i…
Browse files Browse the repository at this point in the history 
…ntersection (#2)

* correct input list algo in Rcpp volume.cpp

* Improvements in rounding and rotating R functions (GeomScale#76)

* add seed in rotate_polytope and update Rd files

* update R rounding function to return the linear map

* fix rounding in V-poly

* overloading rotating() function and remove unused comments

* implement zonotope volume computation with cb algorithm and hpoly in MMC

* fix all bugs in new cooling_hpoly algorithm and improve test function

* all tests for zonotopes with c algortihm passed successfully

* implmentation of diameter computation of all convex bodies except Vpoly-intersection

* add Vpolyintersection diameter computation and volume approximation

* new implmentation of rounding

* add rounding in zonotope volume approximation with Hpoly in MMC

* dispatching in ratio estimation

* update ratio estimation functions

* implmemrnt dispatching for diameter computation

* improve ratio estimation code structure
  • Loading branch information
@TolisChal @vissarion
TolisChal authored and vissarion committed May 18, 2020
1 parent aabd28e commit 2d29e2b
Show file tree
Hide file tree
Showing 24 changed files with 1,245 additions and 605 deletions.
19 changes: 18 additions & 1 deletion include/convex_bodies/ballintersectconvex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,10 @@ class BallIntersectPolytope {
Polytope P;
CBall B;
public:
typedef typename Polytope::MT MT;
typedef typename Polytope::VT VT;
typedef typename CBall::NT NT;
typedef typename CBall::PointType PointType;
typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
NT diameter;

BallIntersectPolytope() {}
Expand Down Expand Up @@ -53,6 +54,18 @@ class BallIntersectPolytope {
return diameter;
}

MT get_mat() const {
return P.get_mat();
}

MT get_T() const {
return P.get_mat();
}

MT get_vec() const {
return P.get_vec();
}

int is_in(PointType &p) const
{
if (B.is_in(p)==-1)
Expand All @@ -77,6 +90,10 @@ class BallIntersectPolytope {
comp_diam(diam);
}

NT radius() {
return B.radius();
}

std::pair<NT,NT> line_intersect(PointType &r, PointType &v) const
{

Expand Down
4 changes: 4 additions & 0 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,10 @@ class HPolytope{
diam = 4.0 * std::sqrt(NT(_d)) * cheb_rad;
}

NT radius() const {
return NT(0);
}

void init(const unsigned int dim, const MT &_A, const VT &_b) {
_d = dim;
A = _A;
Expand Down
67 changes: 47 additions & 20 deletions include/convex_bodies/vpolyintersectvpoly.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,16 +15,18 @@
#include <iterator>
#include <vector>

template <typename VPolytope>
template <typename VPolytope, typename RNGType>
class IntersectionOfVpoly {
public:
typedef typename VPolytope::NT NT;
typedef typename VPolytope::PolytopePoint PolytopePoint;
typedef PolytopePoint Point;
typedef typename VPolytope::PointType PointType;
typedef PointType Point;
typedef typename VPolytope::MT MT;
typedef typename VPolytope::VT VT;
typedef typename VPolytope::rngtype RNGType;
std::vector<Point> vecV;
std::pair<Point,NT> _inner_ball;
NT diameter;
//typedef typename VPolytope::rngtype RNGType;
//std::vector<Point> vecV;
NT rad;
VPolytope P1;
VPolytope P2;
Expand All @@ -36,7 +38,12 @@ class IntersectionOfVpoly {
VPolytope first() { return P1; }
VPolytope second() { return P2; }

int is_in(const Point &p){
std::pair<Point,NT> InnerBall() const
{
return _inner_ball;
}

int is_in(const Point &p) const {
if(P1.is_in(p)==-1)
return P2.is_in(p);
return 0;
Expand Down Expand Up @@ -64,10 +71,22 @@ class IntersectionOfVpoly {
//return 4;
}

std::vector<Point> get_vertices() const {
return vecV;
void set_diameter(const NT &diam) {
diameter = diam;
}

NT get_diameter() const {
return diameter;
}

NT radius() const{
return NT(0);
}

//std::vector<Point> get_vertices() const {
// return vecV;
//}

NT getRad() const {
return rad;
}
Expand All @@ -76,6 +95,10 @@ class IntersectionOfVpoly {
return P1.get_mat();
}

VT get_vec() const {
return P1.get_vec();
}

MT get_T() const {
return P1.get_mat();
}
Expand Down Expand Up @@ -108,8 +131,10 @@ class IntersectionOfVpoly {

bool is_feasible() {
bool empty;
int k = P1.get_mat().rows() + P2.get_mat().rows();
RNGType rng(k);
PointInIntersection<VT>(P1.get_mat(), P2.get_mat(),
get_direction<RNGType, Point, NT>(P1.get_mat().rows() + P2.get_mat().rows()), empty);
GetDirection<Point>::apply(k, rng), empty);
return !empty;
}

Expand All @@ -119,14 +144,15 @@ class IntersectionOfVpoly {
MT V1 = P1.get_mat(), V2 = P2.get_mat();
int k1 = V1.rows(), k2 = V2.rows();
int k = k1 + k2;
RNGType rng(k);
Point direction(k), p(d);
std::vector<Point> vertices;
typename std::vector<Point>::iterator rvert;
bool same;

while(num<d+1){

direction = get_direction<RNGType, Point, NT>(k);
direction = GetDirection<Point>::apply(k, rng);
p = PointInIntersection<VT>(V1, V2, direction, same);

same = false;
Expand All @@ -143,7 +169,8 @@ class IntersectionOfVpoly {

}

return P1.get_center_radius_inscribed_simplex(vertices.begin(), vertices.end());
_inner_ball = P1.get_center_radius_inscribed_simplex(vertices.begin(), vertices.end());
return _inner_ball;

}

Expand Down Expand Up @@ -200,7 +227,7 @@ class IntersectionOfVpoly {

// compute intersection point of ray starting from r and pointing to v
// with the V-polytope
std::pair<NT,NT> line_intersect(const Point &r, const Point &v) {
std::pair<NT,NT> line_intersect(const Point &r, const Point &v) const {

std::pair <NT, NT> P1pair = P1.line_intersect(r, v);
std::pair <NT, NT> P2pair = P2.line_intersect(r, v);
Expand All @@ -212,19 +239,19 @@ class IntersectionOfVpoly {
// compute intersection point of ray starting from r and pointing to v
// with the V-polytope
std::pair<NT,NT> line_intersect(const Point &r, const Point &v, const VT &Ar,
const VT &Av) {
const VT &Av) const {
return line_intersect(r, v);
}


// compute intersection point of ray starting from r and pointing to v
// with the V-polytope
std::pair<NT,NT> line_intersect(const Point &r, const Point &v, const VT &Ar,
const VT &Av, const NT &lambda) {
const VT &Av, const NT &lambda) const {
return line_intersect(r, v);
}

std::pair<NT, int> line_positive_intersect(const Point &r, const Point &v) {
std::pair<NT, int> line_positive_intersect(const Point &r, const Point &v) const {

std::pair<NT, int> P1pair = P1.line_positive_intersect(r, v);
std::pair<NT, int> P2pair = P2.line_positive_intersect(r, v);
Expand All @@ -236,13 +263,13 @@ class IntersectionOfVpoly {
}

std::pair<NT, int> line_positive_intersect(const Point &r, const Point &v, const VT &Ar,
const VT &Av) {
const VT &Av) const {
return line_positive_intersect(r, v);
}


std::pair<NT, int> line_positive_intersect(const Point &r, const Point &v, const VT &Ar,
const VT &Av, const NT &lambda_prev) {
const VT &Av, const NT &lambda_prev) const {
return line_positive_intersect(r, v);//, Ar, Av);
}

Expand All @@ -251,7 +278,7 @@ class IntersectionOfVpoly {
// with the V-polytope
std::pair<NT,NT> line_intersect_coord(const Point &r,
const unsigned int &rand_coord,
const VT &lamdas) {
const VT &lamdas) const {
std::pair <NT, NT> P1pair = P1.line_intersect_coord(r, rand_coord, lamdas);
std::pair <NT, NT> P2pair = P2.line_intersect_coord(r, rand_coord, lamdas);
return std::pair<NT, NT>(std::min(P1pair.first, P2pair.first),
Expand All @@ -265,7 +292,7 @@ class IntersectionOfVpoly {
const Point &r_prev,
const unsigned int &rand_coord,
const unsigned int &rand_coord_prev,
const VT &lamdas) {
const VT &lamdas) const {
return line_intersect_coord(r, rand_coord, lamdas);
}

Expand Down Expand Up @@ -310,7 +337,7 @@ class IntersectionOfVpoly {

void normalize() {}

void compute_reflection (Point &v, const Point &p, const int &facet) {
void compute_reflection (Point &v, const Point &p, const int &facet) const {

if (facet == 1) {
P1.compute_reflection (v, p, facet);
Expand Down
9 changes: 6 additions & 3 deletions include/convex_bodies/vpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,14 +20,14 @@
//min and max values for the Hit and Run functions

// V-Polytope class
template <typename Point, typename RNGType>
template <typename Point>
class VPolytope{
public:
typedef Point PointType;
typedef typename Point::FT NT;
typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
typedef RNGType rngtype;
//typedef RNGType rngtype;

private:
MT V; //matrix V. Each row contains a vertex
Expand Down Expand Up @@ -75,7 +75,6 @@ class VPolytope{
return 0;
}


// compute the number of facets of the cyclic polytope in dimension _d with the same number of vertices
// this is an upper bound for the number of the facets from McMullen's Upper Bound Theorem
unsigned int upper_bound_of_hyperplanes() const {
Expand Down Expand Up @@ -116,6 +115,10 @@ class VPolytope{
return V;
}

NT radius() const {
return NT(0);
}

void init(const unsigned int &dim, const MT &_V, const VT &_b) {
_d = dim;
V = _V;
Expand Down
33 changes: 21 additions & 12 deletions include/convex_bodies/zonoIntersecthpoly.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class ZonoIntersectHPoly {
typedef typename HPolytope::NT NT;
typedef typename HPolytope::VT VT;
typedef typename HPolytope::MT MT;
typedef typename Zonotope::PointType Point;
typedef typename Zonotope::PointType PointType;
NT diameter;

ZonoIntersectHPoly() {}
Expand All @@ -28,7 +28,7 @@ class ZonoIntersectHPoly {
Zonotope first() const { return Z; }
HPolytope second() const { return HP; }

int is_in(const Point &p){
int is_in(const PointType &p) const {
if(HP.is_in(p)==-1)
return Z.is_in(p);
return 0;
Expand All @@ -54,6 +54,10 @@ class ZonoIntersectHPoly {
return HP.get_mat();
}

MT get_T() const {
return Z.get_mat();
}

MT get_vec() const {
return HP.get_vec();
}
Expand All @@ -71,23 +75,28 @@ class ZonoIntersectHPoly {
return diameter;
}

std::pair<NT,NT> line_intersect(Point &r, Point &v) const {
std::pair<PointType,NT> InnerBall() const
{
return Z.InnerBall();
}

std::pair<NT,NT> line_intersect(PointType &r, PointType &v) const {

std::pair <NT, NT> polypair = HP.line_intersect(r, v);
std::pair <NT, NT> zonopair = Z.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, zonopair.first),
std::max(polypair.second, zonopair.second));
}

std::pair<NT,NT> line_intersect(Point &r, Point &v, VT &Ar,
std::pair<NT,NT> line_intersect(PointType &r, PointType &v, VT &Ar,
VT &Av) const {
std::pair <NT, NT> polypair = HP.line_intersect(r, v, Ar, Av);
std::pair <NT, NT> zonopair = Z.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, zonopair.first),
std::max(polypair.second, zonopair.second));
}

std::pair<NT,NT> line_intersect(Point &r, Point &v, VT &Ar,
std::pair<NT,NT> line_intersect(PointType &r, PointType &v, VT &Ar,
VT &Av, NT &lambda_prev) const {
std::pair <NT, NT> polypair = HP.line_intersect(r, v, Ar, Av, lambda_prev);
std::pair <NT, NT> zonopair = Z.line_intersect(r, v);
Expand All @@ -96,7 +105,7 @@ class ZonoIntersectHPoly {
}

//First coordinate ray shooting intersecting convex body
std::pair<NT,NT> line_intersect_coord(Point &r,const unsigned int &rand_coord,
std::pair<NT,NT> line_intersect_coord(PointType &r,const unsigned int &rand_coord,
VT &lamdas) const {

std::pair <NT, NT> polypair = HP.line_intersect_coord(r, rand_coord, lamdas);
Expand All @@ -106,7 +115,7 @@ class ZonoIntersectHPoly {
}


std::pair<NT, int> line_positive_intersect(Point &r, Point &v, VT &Ar,
std::pair<NT, int> line_positive_intersect(PointType &r, PointType &v, VT &Ar,
VT &Av) const {

std::pair <NT, int> polypair = HP.line_positive_intersect(r, v, Ar, Av);
Expand All @@ -118,7 +127,7 @@ class ZonoIntersectHPoly {
return std::pair<NT, int>(std::min(polypair.first, zonopair.first), facet);
}

std::pair<NT, int> line_positive_intersect(Point &r, Point &v, VT &Ar,
std::pair<NT, int> line_positive_intersect(PointType &r, PointType &v, VT &Ar,
VT &Av, NT &lambda_prev) const {
std::pair <NT, int> polypair = HP.line_positive_intersect(r, v, Ar, Av, lambda_prev);
std::pair <NT, int> zonopair = Z.line_positive_intersect(r, v, Ar, Av);
Expand All @@ -130,8 +139,8 @@ class ZonoIntersectHPoly {
}

//Not the first coordinate ray shooting intersecting convex body
std::pair<NT,NT> line_intersect_coord(Point &r,
const Point &r_prev,
std::pair<NT,NT> line_intersect_coord(PointType &r,
const PointType &r_prev,
const unsigned int &rand_coord,
const unsigned int &rand_coord_prev,
VT &lamdas) const {
Expand All @@ -142,14 +151,14 @@ class ZonoIntersectHPoly {
std::max(polypair.second, zonopair.second));
}

std::pair<NT,NT> query_dual(Point &p, const unsigned int &rand_coord) const {
std::pair<NT,NT> query_dual(PointType &p, const unsigned int &rand_coord) const {
std::pair <NT, NT> polypair = HP.query_dual(p, rand_coord);
std::pair <NT, NT> zonopair = Z.line_intersect_coord(p, rand_coord);
return std::pair<NT, NT>(std::min(polypair.first, zonopair.first),
std::max(polypair.second, zonopair.second));
}

void compute_reflection (Point &v, const Point &p, const int &facet) const {
void compute_reflection (PointType &v, const PointType &p, const int &facet) const {

if (facet == (HP.num_of_hyperplanes()+1)) {
Z.compute_reflection(v, p, facet);
Expand Down
Loading

0 comments on commit 2d29e2b

Please sign in to comment.