/* Copyright (c) 1997-2006
Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
http://www.math.tu-berlin.de/polymake, mailto:polymake@math.tu-berlin.de
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version: http://www.gnu.org/licenses/gpl.txt.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
*/
#ifndef _POLYMAKE_BENEATH_BEYOND_H
#define _POLYMAKE_BENEATH_BEYOND_H "$Project: polymake $$Id: beneath_beyond.h 7315 2006-04-02 21:37:53Z gawrilow $"
#include <linalg.h>
#include <Vector.h>
#include <Matrix.h>
#include <ListMatrix.h>
#include <IncidenceMatrix.h>
#include <Graph.h>
#include <Bitset.h>
#include <Set.h>
#include <Array.h>
#include <list>
namespace polymake { namespace polytope {
/** Encapsulating the beneath-beyond algorithm.
One instance is used to solve exactly one convex-hull problem.
@tmplparam E numerical type of the coordinates
*/
template <typename E>
class beneath_beyond_algo {
public:
/** The only constructor, initializing the whole stuff.
@param points_arg rows are points with homogeneous coordinates
@param already_VERTICES_arg true if it is already known that there are no redundant points in the input
*/
beneath_beyond_algo(const Matrix<E>& points_arg, bool already_VERTICES_arg)
: points(points_arg),
already_VERTICES(already_VERTICES_arg), generic_position(already_VERTICES_arg),
facet_normals_valid(false),
AH(unit_matrix<E>(points.cols())),
interior_points(!already_VERTICES ? points.rows() : 0),
free_facet(-1),
vertices_this_step(!already_VERTICES ? points.rows() : 0),
interior_points_this_step(!already_VERTICES ? points.rows() : 0)
{
#ifdef POLY_DEBUG
const char *debug_env=getenv("POLY_DEBUG");
if (!debug_env) {
debug=do_nothing;
} else if (!*debug_env || !strcmp(debug_env, "dump")) {
debug=do_dump;
} else if (!strcmp(debug_env, "full_dump")) {
debug=do_full_dump;
} else if (!strcmp(debug_env, "check")) {
debug=do_check;
} else {
throw std::runtime_error(std::string("beneath_beyond_algo - unknown debug keyword: ") + debug_env);
}
#endif
}
/** The main and only method which performs the computation.
@param perm an end-sensitive input iterator over a sequence of cardinal numbers 0..#points-1
or a permutation thereof. It determines the order in which the points are successively
added, and hence the resulting placing triangulation.
*/
template <typename Iterator>
void compute(Iterator perm);
// The rest public methods must be called after compute(); they merely retrieve the computation results.
// The names are self-explanatory.
Matrix<E> getFacets() const;
IncidenceMatrix<> getVertexFacetIncidence() const;
Matrix<E> getAffineHull() const;
Matrix<E> getVertices() const;
Graph<> getDualGraph() const;
Array< Set<int> > getTriangulation() const;
bool getGenericPosition() const;
protected:
// connects a facet with a triangulation simplex
struct incident_simplex {
const Set<int>* simplex; // an element of triangulation
int opposite_vertex; // the only vertex NOT belonging to the facet
incident_simplex(const Set<int>& simplex_arg, int vertex_arg)
: simplex(&simplex_arg), opposite_vertex(vertex_arg) { }
};
// description of a facet; stored as a node attribute of the dual graph
struct facet_info {
Vector<E> normal; // normal vector, directed inside the polyhedron
E sqr_normal; // sqr(normal)
int orientation; // sign(normal * current-point-to-be-added) during one algo step
// =0 : facet is incident, <0 : facet is violated (and will die)
Set<int> vertices; // ... incident to the facet
// simplices from the polytope triangulation contributing to the triangulation of this facet
typedef std::list<incident_simplex> simplex_list;
simplex_list simplices;
// compute the normal vector etc. assuming the full-dimensional case (Affine Hull is empty)
void coord_full_dim(const beneath_beyond_algo&);
// compute the normal vector etc. in the low-dimensional case
void coord_low_dim(const beneath_beyond_algo&);
// check the intersection of triangulation simplices from [s, s_end) with the facet,
// include those with a (d-1)-face ompletely belonging to the facet in its simplex_list
template <typename Iterator>
void add_incident_simplices(Iterator s, Iterator s_end);
// enables efficient memory management in the dual graph
friend void relocate(facet_info* from, facet_info* to) {
pm::relocate(&from->normal, &to->normal);
pm::relocate(&from->sqr_normal, &to->sqr_normal);
to->orientation=from->orientation;
pm::relocate(&from->vertices, &to->vertices);
pm::relocate(&from->simplices, &to->simplices);
}
};
const Matrix<E>& points;
bool already_VERTICES, generic_position, facet_normals_valid;
// the node attributes are facet descriptions
// the edge attributes are ridges, that is, vertices incident to both facets
typedef Graph<facet_info, Set<int> > DualGraph;
DualGraph dual_graph;
ListMatrix< SparseVector<E> > AH, // affine hull
facet_nullspace; // its affine nullspace - is not computed as long as
// the consumed input vertices comprise a simplex
Bitset interior_points, // indices of points which are not vertices
active_facets;
typedef pm::IndexedSubset<Nodes<DualGraph>&, const Bitset&> ActiveFacets;
typedef std::list< Set<int> > Triangulation;
Triangulation triangulation;
int triang_size; // = triangulation.size();
int valid_facet, // a facet where to start the visibility border search
free_facet; // head of a list of empty (deleted) facets to be recycled
// These are working variables valid within one algo step.
// We define them as instance variables nevertheless to avoid the repeating allocation and deallocation.
Bitset vertices_this_step, // points proved to be non-redundant
interior_points_this_step, // points that could be redundant
visited_facets; // facets seen
// accumulates the non-redundant points; is filled until the polytope turns out to be full-dimensional
Set<int> vertices_so_far;
// get an index of a new facet.
// May force the dual graph to resize, invalidating all references to node attributes!
int new_facet();
// release the most resources in the facet description, put it into recycling list
void delete_facet(int f);
/// create the start configuration: a line segment connecting two points
/// @param p1,p2 row indices in the points matrix
void start_with_points(int p1, int p2);
// add the next point, given by the row index.
// Recalculates the affine hull. If its dimension did not decrease, delegates the rest work to add_point_full_dim()
void add_point_low_dim(int p);
// The first phase of the step: looking for a facet violated by point p. If found, calls update_facets()
void add_point_full_dim(int p);
// helper function for add_point_full_dim
int descend_to_violated_facet(int f, int p);
// helper function
void facet_normals_low_dim();
/// The main phase of the step: detect all facets to be deleted, create new facets and simplices
/// @param f the index of the first facet violated by or incident with point p
void update_facets(int f, int p);
/// The final phase of the step: create new edges in the dual graph
/// @param facets facets incident to the point p
template <typename FacetSubset>
void update_dual_graph(FacetSubset facets);
#ifdef POLY_DEBUG
void dump() const;
void dump_p(int p) const;
void check_f(int f, int last_p) const;
void check_p(int p) const;
void check_fp(int f_index, const facet_info& f, int p, std::ostringstream& errors) const;
enum debug_kind { do_nothing, do_check, do_dump, do_full_dump };
debug_kind debug;
Bitset points_so_far;
#endif
};
template <typename E> inline
Matrix<E> beneath_beyond_algo<E>::getFacets() const
{
return Matrix<E>(dual_graph.nodes(), points.cols(),
entire(attach_member_accessor(nodes(dual_graph), ptr2type<facet_info, Vector<E>, &facet_info::normal>())));
}
template <typename E> inline
Matrix<E> beneath_beyond_algo<E>::getAffineHull() const
{
return AH;
}
template <typename E> inline
Matrix<E> beneath_beyond_algo<E>::getVertices() const
{
return points.minor(~interior_points,All);
}
template <typename E> inline
IncidenceMatrix<> beneath_beyond_algo<E>::getVertexFacetIncidence() const
{
IncidenceMatrix<> VIF(dual_graph.nodes(), points.rows(),
attach_member_accessor(nodes(dual_graph), ptr2type<facet_info, Set<int>, &facet_info::vertices>()).begin());
if (already_VERTICES) return VIF;
return VIF.minor(All,~interior_points);
}
template <typename E> inline
Graph<> beneath_beyond_algo<E>::getDualGraph() const
{
return Graph<>(dual_graph);
}
template <typename E> inline
Array< Set<int> > beneath_beyond_algo<E>::getTriangulation() const
{
return Array< Set<int> >(triang_size, triangulation.rbegin());
}
template <typename E> inline
bool beneath_beyond_algo<E>::getGenericPosition() const
{
return generic_position;
}
} } // end namespace polymake
#endif // _POLYMAKE_BENEATH_BEYOND_H
// Local Variables:
// mode:C++
// c-basic-offset:3
// End:
syntax highlighted by Code2HTML, v. 0.9.1