1 /* RegularTriangulationWithCddlib.cpp -- Regular triangulation using CDDLIB
2
3 Copyright 2006 Matthias Koeppe
4
5 This file is part of LattE.
6
7 LattE is free software; you can redistribute it and/or modify it
8 under the terms of the version 2 of the GNU General Public License
9 as published by the Free Software Foundation.
10
11 LattE is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with LattE; if not, write to the Free Software Foundation,
18 Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
19 */
20
21 #include <cassert>
22 #include <vector>
23 #include "latte_cddlib.h"
24 #include "latte_gmp.h"
25 #include "latte_random.h"
26 #include "RegularTriangulationWithCddlib.h"
27
28 using namespace std;
29
30 listCone *
cone_from_ray_set(vector<listVector * > & rays,set_type ray_set,Vertex * vertex)31 cone_from_ray_set(vector<listVector *> &rays,
32 set_type ray_set,
33 Vertex *vertex)
34 {
35 listCone *c = createListCone();
36 c->vertex = new Vertex(*vertex);
37 vector<listVector *>::iterator i;
38 int j;
39 for (i = rays.begin(), j = 1; i!=rays.end(); ++i, j++) {
40 if (set_member(j, ray_set)) {
41 c->rays = new listVector((*i)->first, c->rays);
42 }
43 }
44 return c;
45 }
46
47 void
triangulate_cone_with_cddlib(listCone * cone,BarvinokParameters * Parameters,height_function_type height_function,void * height_function_data,int cone_dimension,ConeConsumer & consumer)48 triangulate_cone_with_cddlib(listCone *cone,
49 BarvinokParameters *Parameters,
50 height_function_type height_function,
51 void *height_function_data,
52 int cone_dimension,
53 ConeConsumer &consumer)
54 {
55 // Copy rays into an array, so we can index them.
56 vector<listVector *> rays = ray_array(cone);
57 /* Create a matrix from the rays, with 2 extra coordinates in
58 front: The 0th for homogenization, the 1st for lifting. */
59 dd_MatrixPtr matrix
60 = rays_to_cddlib_matrix(cone->rays, Parameters->Number_of_Variables,
61 /* num_homogenization_vars: */ 2,
62 /* num_extra_rows: */ 1);
63 int num_rays = lengthListVector(cone->rays);
64 assert(num_rays + 1 == matrix->rowsize);
65 /* Extra row: `vertical' ray -- This kills all upper facets.
66 See Verdoolaege, Woods, Bruynooghe, Cools (2005). */
67 dd_set_si(matrix->matrix[num_rays][1], 1);
68
69 /* Compute random lifting. */
70 int i;
71 for (i = 0; i<num_rays; i++) {
72 height_function(matrix->matrix[i][1], rays[i]->first, height_function_data);
73 }
74 #if 0
75 /* Output of the file -- for debugging. */
76 {
77 FILE *mf = fopen("lifted_cone_for_triangulation", "w");
78 dd_WriteMatrix(mf, matrix);
79 fclose(mf);
80 }
81 #endif
82 /* Compute facets by double-description method. */
83 dd_ErrorType error;
84 dd_PolyhedraPtr poly = dd_DDMatrix2Poly(matrix, &error);
85 check_cddlib_error(error, "cone_to_cddlib_polyhedron");
86 dd_MatrixPtr inequalities = dd_CopyInequalities(poly);
87 assert(inequalities->representation == dd_Inequality);
88 int num_inequalities = inequalities->rowsize;
89 /* For each computed facet, obtain the set of input rays that are
90 incident with the facet. */
91 dd_SetFamilyPtr incidence = dd_CopyIncidence(poly);
92 assert(incidence->setsize == num_rays + 1);
93 assert(incidence->famsize == num_inequalities);
94 /* Walk through all facets. */
95 for (i = 0; i<num_inequalities; i++) {
96 if (set_member(i + 1, inequalities->linset)) {
97 /* We ignore equalities. */
98 }
99 else {
100 /* If the extra ray is incident, we have a facet that is a wall
101 of the half-infinite cylinder over the cone -- ignore it.
102 Otherwise, the facet gives a simplicial cone of the
103 triangulation. */
104 if (!set_member(num_rays + 1 /* 1-based */,
105 incidence->set[i])) {
106 listCone *c = cone_from_ray_set(rays, incidence->set[i], cone->vertex);
107 /* Is a cone of the triangulation -- check it is simplicial */
108 int c_num_rays = set_card(incidence->set[i]);
109 if (c_num_rays > cone_dimension && !Parameters->nonsimplicial_subdivision) {
110 cerr << "Found non-simplicial cone (" << c_num_rays << "rays) "
111 << "in polyhedral subdivision, triangulating it recursively." << endl;
112 /* In the refinement step, always fall back to using a
113 random height vector. */
114 triangulate_cone_with_cddlib(c, Parameters,
115 random_height, &Parameters->triangulation_max_height,
116 cone_dimension, consumer);
117 freeCone(c);
118 }
119 else if (c_num_rays < cone_dimension) {
120 cerr << "Lower-dimensional cone in polyhedral subdivision, should not happen."
121 << endl;
122 abort();
123 }
124 else {
125 consumer.ConsumeCone(c);
126 }
127 }
128 }
129 }
130 dd_FreeMatrix(inequalities);
131 dd_FreeSetFamily(incidence);
132 dd_FreeMatrix(matrix);
133 dd_FreePolyhedra(poly);
134 }
135
136 void
random_regular_triangulation_with_cddlib(listCone * cone,BarvinokParameters * Parameters,ConeConsumer & consumer)137 random_regular_triangulation_with_cddlib(listCone *cone,
138 BarvinokParameters *Parameters,
139 ConeConsumer &consumer)
140 {
141 triangulate_cone_with_cddlib(cone, Parameters, random_height, &Parameters->triangulation_max_height,
142 Parameters->Number_of_Variables, consumer);
143 }
144
145 void
refined_delone_triangulation_with_cddlib(listCone * cone,BarvinokParameters * Parameters,ConeConsumer & consumer)146 refined_delone_triangulation_with_cddlib(listCone *cone,
147 BarvinokParameters *Parameters,
148 ConeConsumer &consumer)
149 {
150 triangulate_cone_with_cddlib(cone, Parameters, delone_height, NULL,
151 Parameters->Number_of_Variables, consumer);
152 }
153