From 8fc34eb74e4eac62e9e1d9b2b68d67add6bd50db Mon Sep 17 00:00:00 2001
From: Sven Verdoolaege <skimo@kotnet.org>
Date: Fri, 13 Apr 2007 01:17:27 +0200
Subject: [PATCH] volume.c: volume_triangulate: use vertex instead of
 barycenter by default

Proposed by Benoit Meister.
---
 barvinok/options.h |  7 +++++--
 options.c          | 15 ++++++++++-----
 volume.c           | 35 +++++++++++++++++++++++++++++++----
 3 files changed, 46 insertions(+), 11 deletions(-)

diff --git a/barvinok/options.h b/barvinok/options.h
index e891750..fdd82f4 100644
--- a/barvinok/options.h
+++ b/barvinok/options.h
@@ -56,7 +56,10 @@ struct barvinok_options {
     #define	BV_APPROX_SCALE_NARROW2	(1 << 2)
     #define	BV_APPROX_SCALE_CHAMBER	(1 << 3)
     int			    scale_flags;
-    int			    volume_triangulate_lift;
+    #define	BV_VOL_LIFT		0
+    #define	BV_VOL_VERTEX		1
+    #define	BV_VOL_BARYCENTER	2
+    int			    volume_triangulate;
 
     /* basis reduction options */
     #define	BV_GBR_NONE	0
@@ -88,7 +91,7 @@ void barvinok_options_free(struct barvinok_options *options);
 #define BV_OPT_POLAPPROX	261
 #define BV_OPT_APPROX		262
 #define BV_OPT_SCALE		263
-#define BV_OPT_NO_LIFT  	264
+#define BV_OPT_VOL		264
 #define BV_OPT_LAST		264
 
 #define BV_GRP_APPROX		1
diff --git a/options.c b/options.c
index 9fd6b1f..f332a09 100644
--- a/options.c
+++ b/options.c
@@ -72,7 +72,7 @@ struct barvinok_options *barvinok_options_new_with_defaults()
     options->polynomial_approximation = BV_APPROX_SIGN_NONE;
     options->approximation_method = BV_APPROX_NONE;
     options->scale_flags = 0;
-    options->volume_triangulate_lift = 1;
+    options->volume_triangulate = BV_VOL_VERTEX;
 
 #ifdef HAVE_LIBGLPK
     options->gbr_lp_solver = BV_GBR_GLPK;
@@ -118,8 +118,8 @@ static struct argp_option approx_argp_options[] = {
 	"method to use in polynomial approximation [default: drop]" },
     { "scale-options",	    BV_OPT_SCALE,
 	"fast|slow,narrow|narrow2,chamber",	0 },
-    { "no-lift",	    BV_OPT_NO_LIFT,	    NULL,	    0,
-	"don't perform lifting triangulation in volume computation" },
+    { "volume-triangulation",	    BV_OPT_VOL,	    "lift|vertex|barycenter",    0,
+	"type of triangulation to perform in volume computation [default: vertex]" },
     { 0 }
 };
 
@@ -190,8 +190,13 @@ static error_t approx_parse_opt(int key, char *arg, struct argp_state *state)
 		argp_error(state, "unknown suboption '%s'\n", subopt);
 	    }
 	break;
-    case BV_OPT_NO_LIFT:
-	options->volume_triangulate_lift = 0;
+    case BV_OPT_VOL:
+	if (!strcmp(arg, "lift"))
+	    options->volume_triangulate = BV_VOL_LIFT;
+	else if (!strcmp(arg, "vertex"))
+	    options->volume_triangulate = BV_VOL_VERTEX;
+	else if (!strcmp(arg, "barycenter"))
+	    options->volume_triangulate = BV_VOL_BARYCENTER;
 	break;
     case ARGP_KEY_END:
 	if (options->polynomial_approximation == BV_APPROX_SIGN_NONE &&
diff --git a/volume.c b/volume.c
index d07b0dd..7b62079 100644
--- a/volume.c
+++ b/volume.c
@@ -311,6 +311,19 @@ static Matrix *barycenter(Param_Polyhedron *PP, Param_Domain *D)
     return center;
 }
 
+static Matrix *triangulation_vertex(Param_Polyhedron *PP, Param_Domain *D,
+				    Polyhedron *F)
+{
+    Param_Vertices *V;
+
+    FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
+	return V->Vertex;
+    END_FORALL_PVertex_in_ParamPolyhedron;
+
+    assert(0);
+    return NULL;
+}
+
 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
  * If F is a simplex, then the volume is computed of a recursive pyramid
  * over F with the points already in matrix.
@@ -339,16 +352,24 @@ static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
     evalue mone;
     Matrix *center;
     unsigned cut_MaxRays = options->MaxRays;
-    unsigned nparam = D->Domain->Dimension;
+    unsigned nparam = PP->V->Vertex->NbColumns-2;
+    Vector *v = NULL;
 
     POL_UNSET(cut_MaxRays, POL_INTEGER);
 
     value_init(mone.d);
     evalue_set_si(&mone, -1, 1);
 
-    center = barycenter(PP, D);
+    if (options->volume_triangulate == BV_VOL_BARYCENTER)
+	center = barycenter(PP, D);
+    else
+	center = triangulation_vertex(PP, D, F);
     for (j = 0; j < dim; ++j)
 	matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
+    if (options->volume_triangulate == BV_VOL_BARYCENTER)
+	Matrix_Free(center);
+    else
+	v = Vector_Alloc(1+nparam+1);
 
     if (row == 0) {
 	for (j = 0; j < dim; ++j)
@@ -365,6 +386,11 @@ static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 	Param_Domain *FD;
 	if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 	    continue;
+	if (options->volume_triangulate != BV_VOL_BARYCENTER) {
+	    Param_Inner_Product(F->Constraint[j], center, v->p);
+	    if (First_Non_Zero(v->p+1, nparam+1) == -1)
+		continue;
+	}
 	FF = facet(F, j, options->MaxRays);
 	FD = face_vertices(PP, D, F, j);
 	tmp = volume_in_domain(PP, FD, dim, matrix, point,
@@ -380,7 +406,8 @@ static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 	Param_Domain_Free(FD);
     }
 
-    Matrix_Free(center);
+    if (options->volume_triangulate != BV_VOL_BARYCENTER)
+	Vector_Free(v);
 
     for (j = 0; j < dim; ++j) {
 	free_evalue_refs(matrix[row][j]);
@@ -562,7 +589,7 @@ static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
     END_FORALL_PVertex_in_ParamPolyhedron;
 
     if (nbV > (dim-row) + 1) {
-	if (options->volume_triangulate_lift)
+	if (options->volume_triangulate == BV_VOL_LIFT)
 	    vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 					  row, options->MaxRays);
 	else
-- 
2.11.4.GIT